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ABSTRACT 


The  feasibility  of  using  a  three  axis  control,  momentum  bias  system  with  magnetic 
actuators  for  sun-oriented  operations  is  explored.  Relevant  equations  of  motion  are 
developed  for  a  sun-oriented  coordinate  system  and  control  laws  are  developed  for: 
initial  spacecraft  capture  after  launch  vehicle  separation;  reorientation  from 
Earth-oriented  to  a  sun-oriented  operations  mode;  sun-oriented  attitude  control;  and 
momentum  wheel  control. 

Simulations  demonstrating  the  stability  and  time  responsiveness  of  the  system  are 
performed.  Sensor  noise  input  tests  are  performed  to  investigate  the  systems 
susceptibility  to  imperfect  conditions.  Cross  product  of  inertia  effects  are  also  input  to 
test  for  system  instability. 
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I.  INTRODUCTION 


A.  BACKGROUND 

Recently,  a  National  Research  Council  study  group  recommended  the  U.  S. 
government  develop  and  launch  several  satellites  to  monitor  the  sun  and  the  effect  of  solar 
radiation  on  the  Earth's  climate.  The  council's  board  on  global  change  wants  a  series  of 
orbiters  to  provide  long  term,  uninterrupted  and  overlapping  observations  of  the  sun. 
Earth-based  sensors  are  unable  to  detect  precise  fluctuations  in  solar  emissions  because 
the  atmosphere  absorbs  much  of  the  energy.  [Ref.  1] 

A  sun-watching  program  spread  over  many  decades  may  be  precluded  by  a  limited 
access  to  space  due  to  rising  costs  and  limited  launch  capabilities.  Therefore,  the  group 
advocated  the  possible  use  of  civilian  or  military  environmental  satellites  to  help  survey  the 
sun.  Satellites  under  the  control  of  the  National  Oceanic  and  Atmospheric  Administration 
and  the  U.  S.  Air  Force  Defense  Meteorological  Satellite  Program  could  possibly  be 
brought  into  service  for  this  purpose,  but  the  study  states  that  satellites  dedicated  solely 
for  this  program  are  needed  to  avoid  interruption  in  the  data  being  taken.  [Ref  1] 

Candidate  satellites  for  this  program  would  have  to  have  long  endurance,  a 
relatively  high  altitude  to  avoid  any  atmospheric  interference,  and  a  precision  sun-pointing 
accuracy  of  at  least  one  half  degree  for  measuring  instrument  requirements.  A  momentum 
bias  system  with  magnetic  control  actuators  appears  to  be  a  good,  relatively  low  cost 
choice  to  satisfy  the  first  and  last  requirements. 

The  popularity  of  momentum  bias  spacecraft  is  due  primarily  to  their  ability  to 
provide  reasonably  good  attitude  accuracy  (better  than  0.5  degree)  with  relatively 
lightweight  and  inexpensive  attitude  control  hardware.  Gyroscopic  stiffiiess  provided  by 
the  momentum  wheel  is  a  key  feature  of  this  system,  especially  if  the  momentum  bias 
vector  can  be  oriented  toward  or  away  from  the  sun  for  safe  mode  recovery.  The  use  of 
magnetic  attitude  control  for  satellite  applications  has  become  quite  popular  due  to  its 
light  weight,  relatively  low  power  consumption,  and  long  life  expectancy.  Additionally, 
magnetic  attitude  control  systems  provide  a  smoother  torque  output  as  opposed  to  the 
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magnetic  attitude  control  systems  provide  a  smoother  torque  output  as  opposed  to  the 
impulsive  output  of  mass  expulsion  systems  and  also  avoid  catastrophic  failure  due  to  the 
lack  of  moving  parts.  [Ref.  2]  Although  magnetic  control  systems  do  not  produce  as 
much  torque  as  mass  expulsion  or  reaction  wheel  systems,  they  do  provide  a  viable 
option  for  satellite  applications  which  do  not  require  rapid  attitude  reorientations. 

B.  HISTORY 

Initial  use  of  passive  magnetic  control  was  implemented  in  1960  with  the  Transit 
navigational  satellite  program.  The  methodology  was  to  align  one  of  the  satellite  axes 
along  the  direction  of  the  Earth's  magnetic  field.  A  year  later,  the  TRAAC  satellite  was 
launched  with  a  system  that  had  the  capability  of  providing  magnetic  stabilization  on 
command.  [Ref.  3]  Other  spacecraft  which  have  used  magnetic  control  include  the 
Nimbus,  ERTS,  ITOS,  RAE  and  OAO  satellites.  On  the  Nimbus  and  ERTS  spacecraft, 
the  magnets  were  used  to  offset  fixed  imbalances  and  could  be  adjusted  via  ground 
commands.  The  ITOS  system  used  magnetics  for  precession  control,  but  the  system  was 
not  autonomous  and  needed  to  be  controlled  from  the  ground.  In  contrast,  a  closed  loop 
magnetic  system  has  been  used  to  regulate  reaction  wheel  momentum  on  the  OAO 
spacecraft.  [Ref.  2]  Finally,  the  HCCM  satellite  used  a  system  which  is  similar  to  the  one 
proposed  for  use  here  and  the  new  Iridium  system  plans  to  use  a  momentum  bias  with 
magnetic  control. 

C.  PURPOSE 

The  purpose  of  this  thesis  is  to  develop,  and  verify  by  simulation,  the  performance 
of  an  autonomous  attitude  control  system  which  can  serve  the  needs  of  the  sun- watching 
satellite  program.  This  system  should  be  able  to  acquire  the  satellite  after  separation  from 
the  laimch  vehicle,  reorient  the  spacecraft  to  a  sun-pointing  attitude  via  ground  command, 
and  then  perform  accurate  sun-pointing  attitude  control  throughout  the  mission  life.  Most 
attitude  control  applications  that  contained  magnetic  torque  rods  for  control  were  used  in 
an  Earth-oriented  configuration.  This  system  must  be  able  to  perform  in  an 
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-Can  past  magnetic  torque  control  laws  be  modified  to  operate  in  this  new  sun- 
oriented  configuration? 

-How  small  can  Beta  angle  become  before  non-linear  assumptions  break  down 
for  the  system  during  Earth  to  sun  satellite  reorientation?  (Beta  angle  is  defined  in 
Chapter  II,  page  7) 

-Can  the  desired  pointing  accuracies  be  achieved  with  this  system? 

-How  long  will  it  take  to  achieve  initial  capture  and  reorientation  to  the  sun? 


3 


II.  SPACECRAFT  CONFIGURATION  AND  MISSION 

ASSUMPTIONS 


The  spacecraft  bus  configuration  to  be  used  is  extracted  from  the  NPS-STEDI 
project  [Ref.  4]  developed  at  the  Naval  Post  Graduate  School  for  the  Advanced 
Spacecraft  Design  class  (AA  4871).  This  provided  a  baseline  for  equipment  selection, 
disturbance  torque  estimation  and  moments  of  inertia.  The  development  of  this  thesis  is 
not  intended  to  be  an  extension  of  the  NPS-STEDI  project,  but  a  general  exploration  of 
magnetic  attitude  control  for  sun-oriented  operations.  The  NPS-STEDI  spacecraft  uses  a 
momentum  biased  system  to  provide  three-axis  stabilization.  Attitude  information  is 
provided  by  sun  sensors  and  a  single  Earth  sensor  with  magnetometer  data  as  a  coarse 
backup  when  Earth  sensor  data  is  unavailable.  The  selected  Earth  sensor  is  moimted  on 
the  momentum  wheel,  which  provides  a  significant  cost  saving  over  using  separate 
systems.  The  system  to  be  developed  here  will  use  most  of  these  design  criteria,  but  the 
selection  of  Earth  sensor  location  will  be  made  after  system  performance  at  low  Beta 
angles  is  analyzed.  Beta  angle  is  discussed  later. 

The  type  of  configuration  to  be  considered  leaves  the  designer  a  great  deal  of 
options  to  explore.  The  first  decision  should  be  which  direction  to  orient  the 
momentum-bias  vector.  It  could  be  oriented  perpendicular  to  the  solar  ecliptic  which 
means  that  the  momentum  vector  does  not  need  to  be  reoriented  by  one  degree  per  day  to 
follow  the  sun.  In  this  configuration,  wheel  speed  would  be  used  to  point  the  solar  arrays 
and  experiment  toward  the  sun.  Unfortunately,  this  configuration  increases  the  number  of 
earth  sensors  required  for  pitch  angle  error  sensing  as  well  as  the  number  of 
communication  antennae  in  use.  Furthermore,  there  may  be  thermal  cycling  problems 
which  arise  in  this  configuration.  Finally,  the  momentum  vector  would  not  necessarily 
keep  the  solar  arrays  in  a  sim-pointed  direction  during  safe  mode/shut  down.  Another 
approach  would  be  to  orient  the  momentum-bias  vector  toward  or  away  from  the  sun. 

This  gives  the  spacecraft  gyroscopic  stiffiiess  which  would  keep  the  solar  arrays  and  the 
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experiment  aligned  with  the  sun.  This  is  ideal  for  safe  mode/shut  down  operation  as  it 
would  help  to  guarantee  power  to  the  satellite  after  a  malfunction  by  keeping  the  solar 
arrays  pointed  in  the  direction  of  the  sun.  However,  the  momentum  vector  would  have  to 
be  reoriented  by  1  degree  per  day.  Also,  unless  two  Earth  sensors  are  used,  pitch  angle 
sensing  data  may  not  be  available  throughout  the  entire  orbit  period  for  low  Beta  angle 
orbit  configurations.  If  the  Earth  sensor  could  somehow  be  mounted  on  the  sun-facing 
side  of  the  satellite,  then  sensing  data  would  be  available  during  the  critical  period  of 
eclipse  (this  configuration  could  cause  problems  with  the  mounting  of  the  solar  arrays). 
When  pitch  angle  sensing  data  is  lost,  magnetometer  data  could  be  used  to  provide  coarse 
pitch  angle  attitude  sensing  for  control. 

For  this  application,  the  fixed  momentum  bias  vector  is  oriented  along  the  negative 
pitch  axis  and  is  used  to  provide  gyroscopic  stiffness  once  the  spacecraft  is  separated  from 
the  launch  vehicle.  In  addition,  the  wheel's  speed  will  be  used  to  control  pitch  angle.  The 
momentum  wheel  size  is  driven  by  the  need  to  keep  the  momentum  vector  pointed  within 
one  half  degree  of  the  sun  line  throughout  an  assumed  maximum  eclipse  period  of  35 
minutes.  This  maximum  eclipse  time  is  estimated  for  a  low  Beta  angle,  1 :30  am/pm  orbit 
and  will  decrease  as  Beta  angle  increases. 

The  next  thing  to  consider  in  this  type  of  approach,  is  what  type  of  torque  source 
to  use  for  momentum  wheel  desaturation  and  roll/yaw  control.  The  two  types  of  systems 
to  consider  are  thrusters  and  magnetic  torque  rods.  Thrusters  will  not  be  used  because  it 
was  decided  that  no  orbit  control  thruster  system  was  required  for  this  mission,  and  a 
thruster  system  would  increase  weight,  complexity,  space  requirements  and  decrease  the 
satellite's  potential  life-span.  Magnetic  torque  rods  will  be  used  for  momentum  wheel 
desaturation  and  roll/yaw  control  (wheel  speed  will  be  maintained  at  3000  RPM  +/- 10%). 
Roll/yaw  corrections  will  only  be  accomplished  during  sunlight  periods.  Error  signals 
generated  by  sun  sensors  mounted  on  the  negative  Y  face  (sun  pointing)  will  be  fed  back 
to  the  attitude  control  system  to  provide  torque  commands  to  the  magnetic  torque  rods. 
The  torque  rods  will  be  used  to  zero  these  error  signals  over  the  course  of  the  sunlight 
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the  sunlight  periods,  thus  providing  near-zero  spacecraft  pointing  error  at  the  beginning 
of  the  eclipse  period.  Individual  aspects  of  this  design  are  discussed  below  and  additional 
assumptions  are  made. 

A.  ORBIT 

The  orbit  selected  was  a  circular,  sun-synchronous  orbit  with  an  altitude  of  550 
kilometers.  This  leads  to  an  orbit  inclination  angle  of  97.58  degrees  and  an  orbit  period 
of  95.65  minutes.  The  driving  force  for  this  orbit  in  the  NPS-STEDI  project  was  the  need 
for  an  ultraviolet  sunrise  on  every  orbit  so  that  a  complete  atmospheric  density  profile 
could  be  measured  by  the  all  reflection  Michelson  interferometer  [Ref.  4].  To  achieve 
this,  the  initial  longitude  of  the  ascending  node  was  selected  to  be  127.5  degrees  which 
corresponds  to  a  2:30  a.m./p.m.  orbit  when  right  ascension  of  the  sun  in  the  ecliptic  plane 
is  90  degrees  (summer  solstice).  For  this  application,  continuous  sun  coverage  is  desired 
(i.e.  6:00  orbit)  but  orbits  ranging  from  6:00  to  1 :30  will  be  explored  to  test  the  validity  of 
attitude  dynamics  assumptions  at  low  Beta  angles. 

B.  BETA  ANGLE  DEVELOPMENT 

Before  discussing  orbital  reference  frames,  the  topic  of  Beta  angle  and  its 
relationship  to  this  problem  needs  to  be  discussed.  Beta  angle  is  defined  as  the  angle 
between  the  orbit  plane  and  the  sun  line,  which  is  the  line  between  the  center  of  the  Earth 
and  the  center  of  the  sun.  This  angle  can  be  seen  in  Figure  1  and  is  defined  using  Figures 
2  and  3  as 

P  =  arcsin[iS'wS'eCi  +  CuSiSw  —  SuCeSiCw]  (1 ) 

where  S  is  the  sine  of  the  angle,  C  is  the  cosine  of  the  angle,  i  is  the  orbital  inclination 
(97.58  degrees),  e  is  the  Earth  spin  axis  tilt  with  respect  to  the  ecliptic  (23.442  degrees),  u 
is  the  right  ascension  of  the  sun  measured  in  the  ecliptic  plane  and  defined  as 

u=u0+udot*time  (2) 
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with  udot=0. 985648  degrees/mean  solar  day;  and  w  is  the  right  ascension  of  the  line  of 
nodes  defined  as 


w=wO+wdot*time 


(3) 


where  wdot  is  the  nodal  regression  rate  given  by 


wdot=[[-9.9639(R/(R+h))^3/5]Ci]/(l-Ecc^2)^2  (degrees/mean  solar  day)  (4) 


with  Ecc  being  the  eccentricity  (zero  for  a  circular  orbit),  h  being  the  altitude  of  the 
satellite  in  orbit,  i  being  the  orbit  inclination,  and  R  being  the  radius  of  the  Earth  [Ref.  5]. 
The  values  for  uO  and  wO  are  selected  to  vary  Beta  angle  to  test  attitude  dynamics 
assumptions  for  simulation  runs.  The  reason  for  varying  Beta  angle  is  to  generate 
different  Earth  to  sun  rotation  angles  which  will  be  derived  later.  These  rotation  angles 
affect  gravity  gradient  and  aerodynamic  disturbance  torque  inputs  for  the  simulation. 
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degrees 


subsolar 

point 

/ 

✓ 

Equator  v=0 


Figure  1 .  Beta  angle  spatial  relationships  with  Earth  to  Sun  rotation  angles 
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Figure  2.  Napier's  circle  relationships  used  in  Beta  angle  definition 


The  key  to  the  topic  of  a  beta  angle  is  the  fact  that  the  angles  ,  e„ ,  and  Ves  in 
Figure  1  are  dependent  on  beta  angle  and  that  these  angles  are  vital  to  the  coordinate 
transformation  rotations  needed  to  reference  the  Earth-oriented  orbital  frame  to  the 
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Sun-oriented  reference  frame  as  discussed  below.  Using  Figures  1  and  4  along  with 
Napier's  Circle  relationships,  the  Earth  to  sun  angles  can  be  calculated  as 


(j)  ej  =  arcsin[cos  P  cos  ]  (5) 

ee,  =  arctan[^]  (6) 

\\ies  =  arcsin[cos  P  sin  ]  (7) 

where  p  is  defined  as  above  and 

+  (8) 


where  is  a  spatial  factor  necessary  to  orient  the  subsolar  point  in  Figure  1  to  the  true 
anomaly  v  of  the  satellite.  This  is  better  seen  by  saying  that  the  satellite  starts  its  orbit  at 
the  equator  (v=0)  and  proceeds  around  the  Earth  with  angular  rate  cOq.  This  means  that 

v-(x>ot  (9) 

The  quantity  is  also  increasing  at  a  rate  equal  to  cOq,  but  it  starts  =  0)  at  the 
subsolar  point  instead  of  the  equator;  thus  the  need  for  the  offset  v^.  Since  the  value  of  Vs 
is  dependent  on  the  position  of  the  subsolar  point,  it  must  be  a  function  of  the  Earth's 
position  in  its  orbit  about  the  sun.  This  relationship  can  be  developed  using  Figure  3  and 
some  spatial  relationship  variables  p,  q,  and  t.  A  spherical  triangle  showing  a  portion  of 
Figure  3  is  given  in  Figure  5.  Napier's  circle  relationships  are  used  to  solve  for  the  angle 
q  as 


q  =  arccos[-cos(l  80  -  inc)cos  e  +  sin(l  80  -  z>jc)sine  cos  w]  (1 0) 
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Now  using  spherical  triangle  relationships,  p  can  be  solved  for  as 

.  rSin(180-mc)sinw,  /ii\ 

;7  =  arcsin[ - ^ - ]-u  ■  (I  I) 

Using  Equations  (lO)  and  (I  I)  and  the  larger  triangle  of  Figure  5,  t  can  be  solved  for  as 


r  =  arcsin[5i!^]  (12) 

Now  using  the  smaller  triangle  and  Napier's  circle  relationships  of  Figure  5  and  Equation 
(1 2),  Vj  can  be  solved  for  as 

Vs  =  t-  arcsin[tan  P  tan(90  -q)]  ( 1 3) 

This  offset  becomes  important  later,  because  the  Magnetic  field  model  of  the  Earth  is 
affected  if  this  offset  is  not  utilized. 


Figure  4.  Napier's  circle  relationships  for  Earth  to  Sun  rotation  angle  definition 


II 


C.  REFERENCE  FRAMES 

1.  Earth-oriented  Orbital  Reference  Frame 

This  reference  frame  is  the  standard  orbital  coordinates  reference  frame  that  is 
used  in  most  attitude  dynamics  publications.  The  positive  X  axis  of  the  frame  is  oriented 
along  the  satellite's  velocity  vector  Avith  the  positive  Z  axis  located  in  the  orbital  plane  in 
a  nadir  pointing  direction  (i.e.  toward  the  center  of  the  Earth).  The  positive  Y  axis  is 
formed  by  taking  the  cross  product  of  the  Z  and  X  axes.  The  Y  axis  is  perpendicular  to 
the  orbital  plane  and  along  the  negative  orbit  normal. 

2,  Sun-oriented  Orbital  Reference  Frame 

This  reference  frame  is  specially  designed  for  this  orbit  application.  In  this 
reference  frame,  the  negative  Y  axis  is  oriented  in  the  direction  of  the  sun  and  the  plane 
formed  by  the  Y  and  Z  axes  cuts  through  the  center  of  the  Earth.  To  orient  the  spacecraft 
in  this  reference  frame,  two  axis  rotations  are  performed  from  the  Earth-oriented  orbital 
reference  frame  using  the  Earth  to  Sun  rotation  angles  discussed  above.  The  first  rotation 
is  performed  about  the  3-axis  or  Zq  axis  of  Figure  1  using  the  angle  \]!^.  This  forms  a  new 
set  of  axes  called  xl,  x2,  and  x3.  It  can  be  shown  that  the  xl-zl  plane  and  the  yl-zl 
plane  still  cut  through  the  center  of  the  Earth  since  the  zl  axis  is  still  oriented  toward  the 
center  of  the  Earth.  Now  a  rotation  about  the  xl  axis  is  performed  using  the  angle 
This  rotation  serves  to  now  orient  the  negative  Y  axis  in  the  direction  of  the  Sun  as  in 
Figure  1.  Now,  since  the  xl  axis  is  perpendicular  to  the  yl-zl  plane  by  definition,  and 
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since  the  yl-zl  plane  cuts  through  the  center  of  the  Earth  before  the  rotation;  it  must  still 
be  cutting  through  the  center  of  the  Earth,  which  is  the  desired  configuration. 

3.  Body-fixed  Reference  Frame 

This  reference  frame  is  obtained  by  performing  a  standard  3-1-2  (Yaw-Roll-Pitch) 
Euler  axis  rotation  in  the  Sun-oriented  orbital  reference  frame.  The  rotations  about  all 
three  axes  will  be  small  angles,  so  small  angle  approximations  can  be  used  in  the 
definition  of  the  equations  of  motion. 

D.  SPACECRAFT  BUS  DESCRIPTION 

The  spacecraft  bus  is  shown  in  Figure  6.  The  main  body  of  the  bus  has 
dimensions  in  inches  of  19  x  19  x  30.  There  are  three  identical  solar  panels  with  the 
center  panel  mounted  directly  to  the  bus  with  the  other  two  panels  hinged  to  either  side. 
All  sun-monitoring  equipment  will  be  mounted  on  the  top  of  the  satellite.  Figure  6  shows 
a  rectangular  box  which  will  act  as  the  equipment  in  this  case  for  disturbance  torque 
analysis.  The  momentum  wheel  and  co-mounted  Earth  sensor  are  located  at  the  back  of 
the  satellite  toward  the  bottom.  The  Earth  sensor  assembly  protrudes  out  the  back  of  the 
satellite  slightly. 

Inertia  and  center  of  mass  properties  for  the  satellite  were  extracted  from 
Reference  4  and  are  listed  below.  The  reference  point  (0,0,0)  for  all  calculations  is 
located  at  the  center  of  the  base  of  the  satellite  (i.e.  on  the  center  of  the  positive  Z  face  at 
the  bottom  of  the  satellite). 

=  8.363  kg-m^ 
lyy  =  9.284  kg-m^ 

=  4.012  kg-m^  (14) 

=  -0.111  kg-m^ 

I„  =  0.013  kg-m^ 
ly,  =  0.048  kg-m^ 
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COM,  =  .000254  m 


COMy  = -.005842  m  (15) 

COM,  = -.453136  m 


E.  ATTITUDE  DETERMINATION  AND  CONTROL  COMPONENTS 

Table  1  gives  a  listing  of  control  system  components  selected  for  the  design. 

Each  component  will  be  explained  thoroughly  below. 

1.  Sensors 

a.  Earth  Sensor 

Knowing  the  spacecraft  attitude  with  respect  to  the  Earth  is  of  great 
importance  in  order  to  maintain  constant  communications  with  ground  stations.  Since  the 
Earth  appears  to  be  rather  large  for  near-Earth  orbiting  satellites,  most  sensors  are 
designed  to  detect  the  Earth's  horizon. 
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component 

size 

mass  (kg) 

power  (watts) 

comments 

Momentum 

wheel 

9.75  N-m-sec 

5.8 

5 

3000  RPM 
+/- 10% 

Magnetic 

Torque  Rods 

15  A-m^ 

0.45  (x3) 

1.43  (x3) 

dual  coil 

Magnetometers 

+/-0.6  Gauss 

0.25  (x2) 

0.04  (x2) 

2-axis 

Earth  Sensor 

1.9 

3.3 

on  momentum 
wheel  initially 

Sun  Sensors 

1.8 

0.5 

11.35 

13.2 

Table  1.  ADCS  components 


The  smallest  amount  of  horizon  sensing  ambiguity  is  achieved  by  using 
the  spectral  region  in  the  infrared  near  the  15pm  wavelength.  Infrared  detectors  are  not 
affected  by  night  or  the  presence  of  a  terminator.  Furthermore,  they  are  less  susceptible 
to  sunlight  reflections  off  the  surface  of  the  spacecraft.  However,  visible  light  sensors 
have  some  advantages  over  infrared  sensors  including:  lower  cost,  faster  response  time 
and  higher  signal  to  noise  ratios  due  to  increased  visible  radiated  intensity.  [Ref.  6] 

The  main  purpose  of  the  Earth  sensor  for  this  satellite  is  to  provide  pitch 
error  data  which  will  be  used  to  control  pitch  angle  and  keep  the  communications  antenna 
(located  on  the  Earth-facing  side)  nominally  Earth-pointed.  This  will  be  accomplished  by 
keeping  the  pitch/yaw  plane  in  such  an  orientation  that  it  cuts  a  great  circle  arc  on  the 
Earth.  This  guarantees  Earth  sensor  data  to  the  maximum  extent  possible.  Due  to  the 
orbit  configuration  and  attitude  control  approach  selected.  Earth  sensor  data  will  not 
always  be  available.  Using  Figure  7,  one  can  see  the  relationships  between  the  Earth 
sensor  scan  line,  which  is  canted  at  a  45  degree  angle  up  from  the  local  vertical,  and  the 
surface  edge  of  the  Earth.  It  is  easy  to  determine  that  if  the  Earth  to  Sun  roll  rotation 
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angle  ((jj^j)  is  greater  than  22  degrees  in  the  negative  direction,  then  the  scan  line  of  the 
sensor  will  not  hit  the  Earth  and  sensing  data  will  be  unavailable.  Current  simulation 
results  show  that  Earth  sensor  data  will  be  unavailable  for  an  increasing  percentage  of  the 
orbit  as  Beta  angle  decreases.  When  Earth  sensor  data  is  unavailable,  the  magnetometers 
and  an  orbital  elements  model  with  a  magnetic  field  model  will  be  used  to  detect  coarse 
pitch  errors  (magnetometer  sensing  error  is  approximately  5  degrees). 

The  Earth  sensor  is  mounted  on  the  momentum  wheel  and  protrudes  out 
the  positive  pitch  axis  surface  of  the  spacecraft  (anti-sun  face).  Consideration  was  given 
to  mounting  the  Earth  sensor  on  the  sun-facing  side  of  the  spacecraft  for  increased 
sensing  data  during  eclipse  when  roll/yaw  error  signals  are  unavailable  fi-om  the  sun 
sensors,  but  solar  array  development  considerations  and  the  cost  of  an  Earth  sensor 
mounted  separately  from  the  momentum  wheel  can  make  this  option  unattractive.  The 
addition  of  this  extra  Earth  sensor  may  become  necessary  to  provide  sensing  data  in 
eclipse  for  low  Beta  angle  applications,  which  will  be  explored  in  the  simulations. 
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b.  Sun  Sensors 

Sun  sensors  have  been  used  on  nearly  all  spacecraft  which  have  flown. 
Since  the  distance  to  the  sun  is  rather  large  for  Earth  orbiting  satellites,  the  angular  radius 
of  the  sun  is  nearly  independent  of  the  type  of  orbit  and  for  most  applications,  the  sun  can 
be  treated  as  a  point-source.  Furthermore,  the  sun  is  bright  enough  to  allow  the  use  of 
simple  equipment  to  discriminate  among  radiation  sources  vdth  small  power 
requirements.  Generally,  sun  sensors  are  used  for  detecting  sun  presence,  directing  the 
solar  arrays  toward  the  sun  and  providing  the  relative  sun  vector  for  protection  of 
sensitive  equipment  such  as  star  trackers  and  other  sensors.  Sun  sensors  can  have  very 
narrow  to  very  wide  fields  of  view  and  resolutions  from  several  degrees  to  less  than  one 
arc-second.  The  three  basic  types  of  sun  sensors  include:  Analog  sensors  which  provide 
a  continuous  output  which  is  a  ftmction  of  sun  angle;  sun  presence  sensors  which  have  a 
constant  output  when  the  sim  is  present;  and  digital  sensors  which  provide  a  discrete 
output  as  a  function  of  sun  angle  [Ref.  6]. 

High  accuracy  sun  sensors  will  be  used  on  the  sun-facing  side  of  the 
spacecraft  to  provide  sun  line  roll/yaw  error  signals  to  the  attitude  control  system  for 
magnetic  torque  rod  control.  Coarse  sim  sensors  will  also  be  used  on  each  face  of  the 
spacecraft  in  case  a  problem  develops  and  the  spacecraft's  attitude  with  respect  to  the  sun 
needs  to  be  known  for  solar  array  pointing.  These  coarse  sun  sensors  will  also  be  needed 
to  provide  pointing  control  while  reorienting  the  spacecraft  fi’om  its  Earth-oriented 
configuration  during  checkout  to  the  sun  pointing  mode. 

c.  Magnetometers 

Magnetometers  are  very  popular  attitude  determination  sensors  and  have 
been  used  widely  on  spacecraft.  They  are  vector  sensors  because  they  provide  both 
direction  and  magnitude  of  the  ambient  magnetic  field.  Furthermore,  they  are  reliable, 
lightweight,  and  have  low  power  requirements;  they  have  no  moving  parts;  and  they  can 
operate  over  a  wide  temperature  range.  Unfortunately,  magnetometers  are  not  very 
accurate  inertial  attitude  sensors  mainly  due  to  the  lack  of  an  accurate  model  of  the 
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Earth's  magnetic  field.  Also,  the  Earth's  magnetic  field  strength  decreases  as  a  function 
of  the  inverse  cube  of  the  satellite  distance  from  the  center  of  the  Earth.  This  could  result 
in  residual  spacecraft  magnetic  effects  dominating  the  field  strength,  generally  limiting 
the  use  of  magnetometers  to  spacecraft  below  1000  kilometers  (note:  missions  have 
successfully  used  magnetometers  above  this  altitude)  [Ref  6]. 

The  magnetometer  consists  of  a  magnetic  sensor  and  an  electronics  unit 
for  transforming  sensor  measurements  into  a  usable  format.  There  are  two  main  types  of 
magnetometers:  quantum  magnetometers,  which  utilize  atomic  properties;  and  induction 
magnetometers,  which  are  based  on  Faraday's  Law  of  Magnetic  Inductance.  The  latter 
group  is  more  popular  and  will  be  discussed  further  here.  [Ref.  6] 

Faraday's  law  states  that  an  electromotive  force  (EMF)  is  induced  in  a 
conducting  coil  when  it  is  placed  in  a  time-varying  magnetic  flux.  There  are  two  types  of 
induction  magnetometers  called  search-coil  and  fluxgate  magnetometers.  Search-coil 
magnetometers  contain  a  solenoid  coil  of  N  turns  surrounded  by  a  ferromagnetic  core 
with  magnetic  permeability  p,  and  cross-sectional  area  A.  The  EMF  induced  in  the  coil 
when  placed  in  a  magnetic  field  produces  a  voltage  as  a  function  of  the  time  varying 
component  of  the  magnetic  field  which  is  parallel  to  the  solenoid.  This  type  of 
magnetometer  is  mainly  used  on  spin-stabilized  spacecraft  to  provide  precise  phase 
information.  The  fluxgate  magnetometer  uses  two  coils  and  two  saturable  cores  to  detect 
the  presence  of  any  ambient  magnetic  field  by  measuring  harmonic  frequencies  formed 
by  the  time-varying  field. 

A  pair  of  two-axis  magnetometers  will  be  used  to  detect  the  ambient 
magnetic  field  the  spacecraft  is  encountering.  The  data  will  be  used  primarily  to  calculate 
electrical  currents  needed  to  generate  desired  control  torques  with  the  magnetic  torque 
rods  and  also  in  a  backup  roll  to  detect  course  attitude  errors  as  discussed  in  the  Earth 
sensor  section  above.  The  magnetometers  will  be  mounted  in  such  a  way  that  each 
element  detects  a  component  of  the  magnetic  field  along  each  major  body  axis.  These 
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signals  will  then  be  averaged  to  determine  the  magnetic  field  strength  eind  direction.  This 
configuration  will  provide  some  backup  in  case  of  failure  in  one  of  the  components. 

2.  Actuators 

a.  Momentum  Wheel 

Angular  momentum  storage  devices  are  used  on  spacecraft  for  several 
purposes  including;  the  addition  of  stability  to  counter  disturbance  torques,  to  provide  a 
variable  momentum  to  allow  operation  at  one  revolution  per  orbit  for  Earth  oriented 
missions,  to  absorb  cyclic  torques,  and  to  transfer  momentum  to  the  spacecraft  for 
slewing  maneuvers.  The  principle  of  momentum  for  a  spinning  wheel  is  derived  from  the 
relationship 


h  =  /(o 


(16) 


where  /  is  the  moment  of  inertia  of  the  wheel  and  ©is  the  rotational  velocity. 

Momentum  wheels  are  designed  to  provide  a  variable  momentum  'bias'  or 

nonzero  momentum  which  can  be  used  for  storage  capability  about  the  rotational  axis  of 
the  wheel  which  is  generally  fixed  with  respect  to  the  spacecraft.  The  same  momentum 
can  be  achieved  with  a  large  inertia,  slowly  rotating  wheel  and  a  small  inertia,  rapidly 
rotating  wheel.  Design  tradeoffs  usually  favor  the  smaller  wheel  for  the  weight  savings, 
but  the  increased  wheel  speed  can  lead  to  greater  wear  on  the  bearings  which  may  shorten 
the  wheel  lifetime. 

Torque  motors  using  either  direct  or  alternating  current  (DC  or  AC)  are 
used  to  transfer  momentum  between  the  wheel  and  the  spacecraft  body.  The  DC  motors 
are  generally  more  efficient  and  provide  high  torques  at  low  speed,  but  they  do  require 
the  use  of  brushes  which  can  wear  out.  AC  motors  do  not  have  this  drawback,  but  they 
do  require  gearing  systems  due  to  their  high  speed.  Also,  they  are  less  efficient  and 
provide  a  lower  torque.  Tachometers  are  used  to  measure  wheel  speed  for  attitude  control 
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The  main  purpose  of  the  momentum  wheel  on  this  spacecraft  is  to  provide 
gyroscopic  stiffness  along  the  pitch  axis  once  the  spacecraft  is  operational.  It  also  serves 
to  control  pitch  angle  by  changing  wheel  speed  to  generate  the  necessary  control  torques. 
The  amount  of  time  that  the  momentum  bias  vector  needs  to  keep  pointing  accuracy  is 
assumed  to  be  35  minutes;  the  maximum  eclipse  time.  Momentum  wheel  size  is  listed  in 
Table  2  for  several  pointing  accuracies  over  this  maximum  eclipse  period  assuming  worst 
case  disturbance  torques  are  present  throughout  the  eclipse.  This  assumption  provides  for 
some  built  in  safety  margin  for  the  sizing. 


accuracy  after  eclipse  period  (degrees) 

momentum  wheel  size  (N-m-sec) 

0.1 

25 

0.2 

12.5 

0.5 

4.8 

1 

2.4 

2 

1.2 

Table  2.  Momentum  wheel  size  versus  gyroscopic  stiffness  pointing  accuracy 


The  amount  of  torque  output  required  is  driven  by  lyy  and  pitch  angular 
acceleration  (0),  which  is  a  function  of  Beta  angle.  Equations  for  0  are  developed  in 
Reference  4  and  are  listed  here  as 


0 


max 


coo^sin^ 

(1-CTCOS^)^ 


(17) 


where 


4cos^psin  P 
(1+sin^P)^ 


(18) 
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(19) 


1-sin^P 

a  = - :r 

l+sin^P 

5  =  arccos[J^ii|^]  '  (20) 

The  maximum  momentum  change  required  from  the  wheel  is  given  by 

•  •  1  .  cos^  P 

wheel  —  (0  max  “  6  niin)^XV  ~  (siiip  ~  sin  P)Q)  olyy  —  ^jjjp  ®  olyy  (21) 

The  maximum  torque  required  from  the  wheel  is  then  determined  from 

Tmeoi  —  Jyy^max.  (22) 

A  wheel  with  a  momentum  of  9.75  N-m-sec  at  3000  rpm  with  a  torque  rating  of  0.025 
N-m  was  selected.  This  provides  a  pointing  deviation  slightly  less  than  0.5  degrees  over 
the  maximum  eclipse  period  assuming  worst  case  disturbance  torques  are  acting  on  the 
spacecraft  throughout  this  period.  Initial  simulation  results  indicate  that  the  maximum 
pitch  angular  acceleration  can  be  shown  to  be  approximately  2x10'®  rad/secl  This 
means  there  can  be  an  extremely  large  lyy  before  wheel  torque  becomes  a  concern. 
b.  Magnetic  Torque  Rods 

Magnetic  coils  or  electromagnets  are  used  to  generate  torques  necessary  to 
control  spacecraft  attitude  and  momentum.  This  could  include  use  for  overcoming 
disturbance  torques  and  compensating  for  residual  magnetic  effects.  The  magnetic 
moment  generated  by  a  loop  of  wire  possessing  N  turns  is  given  by 

m  =  NIAn  (23) 

where  N  is  the  number  of  turns,  I  is  the  current  through  the  loop,  A  is  the  area  enclosed 
by  the  loop  and  n  is  a  imit  vector  normal  to  the  plane  of  the  loop.  The  positive  sense  of 


21 


the  direction  of  the  magnetic  moment  is  determined  by  the  right  hand  rule.  The  fingers  of 
the  right  hand  are  curled  around  the  turns  of  wire  in  the  direction  of  the  current  and  then 
the  right  thumb  points  in  the  direction  of  the  magnetic  moment. 

The  magnetic  dipole  moment  depends  on  what  type  of  material  is  enclosed 
within  the  turns  of  wire  and  is  expressed  by 


M  =  pm  (24) 

where  p  is  the  permeability  of  the  core  material.  The  permeability  of  free  space  p^  has 
the  value  An  x  10'’N/A^  Thus,  it  is  clear  that  the  values  of  core  permeability,  coil 
configuration,  N  and  A,  and  the  current  values  must  be  selected  appropriately  for  the 
given  application.  Magnets  on  most  satellites  have  used  an  "air"  core.  To  generate 
magnetic  control  torques  the  following  relationship  is  used 

Tma^  =  MxB  (25) 

where  M  is  the  magnetic  dipole  moment  and  B  is  the  magnetic  field  expressed  in  body 
coordinates. 

The  magnetic  torque  rods  will  be  used  for  momentum  wheel  desaturation 
and  roll/yaw  repositioning  to  zero  out  sun  sensor  error  signals.  The  change  in  momentum 
required  to  reorient  the  momentum  vector  by  1  degree  per  day  is  approximately 
0.122  N-m-s.  This  relates  to  a  constant  torque  of  1.41  x  10"^  N-m  throughout  the  day  (or 
2.82  X  lO"*  N-m  on  a  50%  duty  cycle)  .  This  value  is  much  smaller  than  the  maximum 
disturbance  torque,  so  the  latter  will  be  used  to  size  the  magnetic  torque  rods.  Assuming 
a  worst  case  magnetic  field  strength  of  1.5  x  10'^  tesla,  the  required  torque  rod  moment  is 
approximately  2.45  A-m^ .  Magnetometer  sensing  data  will  be  used  to  provide  control 
signals  for  electrical  current  generation  which  control  the  direction  and  amount  of  torque 
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output.  The  selected  torque  rods  have  two  coils,  which  provides  redundancy  in  the 
system,  and  have  a  magnetic  moment  of  15  A-m^. 
c.  Thrusters 

Since  no  orbit  control  is  considered  to  be  necessary,  a  thruster  system 
was  omitted  due  to  a  great  increase  in  weight,  complexity  and  cost. 

F.  SYSTEM  PERFORMANCE 

Once  the  injection  vehicle  has  attained  the  proper  orbit,  the  momentum  wheel  will 
be  spim  up  with  the  injection  vehicle  providing  torque  to  counter  the  increase  in 
momentum.  It  will  take  approximately  7  minutes  to  spin  up  the  wheel  if  injection  vehicle 
control  is  available.  In  the  event  the  launch  vehicle  cannot  provide  stabilization,  the 
wheel  will  be  spun  up  in  accordance  with  the  capture  control  section  in  Chapter  IV. 

When  the  wheel  has  achieved  speed,  the  separation  system  will  activate  (assuming 
injection  vehicle  control)  and  separate  the  satellite  from  the  injection  vehicle.  Once  clear, 
the  solar  arrays  will  deploy  and  operations  will  begin  to  establish  a  stable  nadir-pointing 
attitude.  Worst  case  tip  off  rate  in  pitch  is  expected  to  be  approximately  5.5  degrees  per 
second  and  worst  case  nutation  angle  (due  to  roll  tip-off  torque)  is  expected  to  be  5 
degrees.  The  momentum  wheel  is  designed  with  sufficient  momentum  storage  to  absorb 
the  former  and  the  magnetic  torque  rods  will  be  used  to  damp  out  the  latter. 

When  the  satellite  is  stable  in  a  nadir-pointing  configuration,  system  checks  will 
be  performed.  Next,  ground  commands  will  be  used  to  initiate  the  reorientation  of  the 
momentum  bias  vector  toward  the  sun  with  torque  rod  control.  Following  this,  the 
attitude  control  system  will  maintain  pointing  accuracy  of  0.5  degree  in  roll  and  yaw  and 
5.0  degrees  in  pitch. 

The  allowable  speed  region  for  the  wheel  will  be  3000+/-300  RPM.  Wheel 
desaturation  will  be  performed  as  needed  using  the  magnetic  torque  rods.  Initial 
simulation  results  show  that  wheel  speed  takes  several  orbits  to  reach  desaturation  limits 
and  that  the  magnetic  control  law  is  easily  capable  of  accomplishing  the  desaturation. 
Other  simulation  results  from  MATLAB/SIMULINK  also  indicate  that  the  momentum 
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bias  vector  is  capable  of  limiting  roll/yaw  drift  to  .2  degrees  throughout  the  maximum 
eclipse  period  (35  minutes)  with  no  control  inputs.  If  the  proper  alignment  is  achieved 
before  eclipse  entry,  then  the  spacecraft  should  be  oriented  at  visible  sunrise  so  that  the 
stm  is  within  the  field  of  view  of  the  sun-monitoring  instruments. 

G.  SIMULATION 

An  old  SIMULINK  program,  developed  by  Professor  Leonard  at  the  Naval  Post 
Graduate  School,  Monterey,  California  has  been  upgraded  and  tailored  for  this  particular 
application.  It  includes  simulation  blocks  for  each  disturbance  torque  (solar  pressure, 
gravity  gradient,  and  aerodynamic)  and  a  simple  simulation  for  the  Earth  sensor  and  sun 
sensors.  From  these  inputs,  spacecraft  torques  will  be  calculated  and  the  spacecraft's 
dynamic  properties  will  be  simulated.  A  separate  MATLAB  program  will  be  used  to 
control  the  SIMULINK  program  and  create  output  displays.  SIMULINK  and  MATLAB 
block  diagrams  and  control  codes  are  listed  in  the  appendix. 
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III.  SPACE  ENVIRONMENT  AND  ITS  AFFECT  ON  THE 

SPACECRAFT 

A.  EARTH'S  MAGNETIC  FIELD  DEVELOPMENT 

The  derivation  of  the  Earth's  magnetic  field  is  very  complex  for  both  magnitude 
and  direction.  However,  the  field  can  be  approximated  by  a  large,  tilted  dipole  located  at 
the  Earth's  center  with  the  'south'  end  of  the  dipole  located  in  the  northern  hemisphere. 
The  location  of  the  magnetic  poles  is  constantly  drifting  at  approximately  .014  degrees 
per  year  [Ref.  6].  More  accurate  models  of  the  field  are  available  but  they  require  vast 
computing  resources  to  simulate;  and  they  are  not  required  for  use  in  this  feedback 
control  application.  Thus,  the  dipole  approximation  will  be  used. 

The  plane  which  is  perpendicular  to  the  dipole  is  called  the  magnetic  equator  and 
this  is  where  the  field  is  weakest  in  magnitude.  Figure  8  shows  how  the  magnetic  field 
strength  varies  with  altitude  at  the  magnetic  equator.  It  can  be  shown  that  the  strength  of 
the  dipole  field  decreases  as  a  fimction  of  the  inverse  cube  of  the  distance  from  the  dipole 
center  (i.e.  the  center  of  the  Earth)  [Ref  6]. 


Figure  8.  Magnetic  Field  Intensity  at  Equator  in  (nano-Tesla)  [Ref  6] 
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As  magnetic  latitude  increases  from  zero  to  90  degrees,  the  field  strength  at  a  given 
altitude  doubles  as  shown  in  Figure  9. 


Figure  9.  Relative  Intensity  of  Magnetic  Field  as  a  Function  of  Latitude  [Ref  6] 

For  magnetic  control,  the  Earth's  magnetic  field  needs  to  be  expressed  in 
coordinates  relative  to  the  three  primary  body  axes  of  the  spacecraft.  The  magnetic  field 
components  expressed  in  the  Earth  oriented  orbital  reference  frame  are:  [Ref.  7] 


Bxo  =  ^{CzSiCa  -  SolS^iSe  -  CiCaC\iSE] 

(26) 

Byo  =  ^[-CeCi-C^iSiSE] 

(27) 

Bzo  =  ^[ICzSiSa  +  2Ca5p58  -  2a5aCp5'e] 

(28) 

Where  C  is  the  cosine  of  the  angle,  S  is  the  sine  of  the  angle,  i  is  the  orbit  inclination,  R  is 
the  radius  of  the  satellite  orbit,  Mis  the  Earth's  Dipole  strength  (7.943  x  10'^  NmVAmp), 
s  is  the  dipole  tilt  with  respect  to  the  Earth's  spin  axis  (1 1 .7  degrees),  p  is  the  Earth's 
rotation  with  respect  to  the  ascending  node  (varies  360  degrees  in  one  day),  and  a  is  the 
orbit  angle  with  respect  to  the  ascending  node  (i.e.  the  angular  position  of  the  satellite 
measured  counterclockwise  around  the  orbit  from  the  ascending  node). 

To  express  these  magnetic  field  components  (Bo*,  Boy,  and  Bq^)  in  the  spacecraft 
body  axis  reference  frame,  two  coordinate  transformations  are  needed.  The  first 
transforms  the  magnetic  field  from  the  Earth-oriented  orbital  reference  frame  to  the 
sun-oriented  orbital  reference  frame,  while  the  second  is  a  standard  3-1-2  Euler  rotation 
which  leads  to  the  body-fixed  reference  frame  orientation.  Once  the  field  components  are 
represented  in  body-fixed  coordinates,  they  can  be  used  in  the  equations  of  motion  for 
both  disturbance  torque  analysis  and  magnetic  control  torque  inputs. 

B.  TORQUES  WHICH  AFFECT  THE  SPACECRAFT 

To  be  able  to  analyze  the  motion  of  the  spacecraft  on  orbit  one  must  be  able  to 
model  all  torques  which  act  upon  the  spacecraft  whether  they  are  generated  by  the 
environment  or  by  the  spacecraft  attitude  control  system.  These  torques  can  either  be 
categorized  as  cyclic,  varying  sinusoidally  over  an  orbit,  or  secular,  accumulating  with 
time  and  not  averaged  over  an  orbit.  Numerical  integration  of  Euler's  equations  requires 
that  all  torques  be  modeled  as  a  function  of  time  and  the  spacecraft's  attitude.  For  this 
analysis,  the  major  environmental  disturbance  torques  to  be  considered  are:  aerodynamic 
drag,  solar  pressure,  residual  magnetic  dipole  and  Earth's  gravitational.  These 
disturbance  torques  would  quickly  reorient  the  spacecraft  if  not  resisted  in  some  way  by 
the  attitude  control  system.  The  active  control  torques  to  be  used  are  generated  by  the 
momentum  wheel  and  magnetic  coils. 

1.  Disturbance  Torques 

Before  sizing  the  ADCS  components,  one  must  know  the  maximum  possible 
disturbance  torques  that  the  spacecraft  might  experience  throughout  its  lifetime.  As 
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discussed  above,  the  major  categories  of  disturbance  torques  are:  aerodynamic,  solar 
pressure,  residual  magnetic  dipole  and  Earth's  gravitational.  The  type  of  torque  which  is 
the  largest  for  this  type  of  low  altitude  application  is  the  aerodynamic  torque.  Also,  since 
there  are  several  techniques  available  to  limit  the  effects  of  residual  magnetic  dipole 
disturbances,  this  category  will  be  ignored  for  the  purpose  of  simplicity  in  the  simulation. 
Using  References  6,  8  and  9,  each  disturbance  torque  is  analyzed  to  determine  its 
magnitude  and  which  axes  of  the  spacecraft  it  affects, 
a.  Aerodynamic 

For  lower  altitude  orbits,  the  interaction  of  the  spacecraft  with  the  Earth's 
upper  atmosphere  causes  a  disturbance  torque  about  the  spacecraft  center  of  mass.  The 
force  due  to  the  impact  of  atmospheric  molecules  on  the  spacecraft  surface  can  be 
modeled  as  an  elastic  impact  without  reflection.  The  incident  particle's  energy  is 
generally  completely  absorbed  [Ref.  6].  The  incremental  force  dF  on  a  surface  element 
area  dA  with  outward  normal  N  is  given  by 

t/F  =  -|Cz)pF^(N-V)Vfi4  (29) 

where  V  is  the  vmit  vector  along  the  translational  velocity  of  the  individual  surface 

element  relative  to  the  incident  stream  of  particles,  Cq  is  the  coefficient  of  drag  (assumed 
to  be  2.0  if  no  measured  value  is  available)  and  p  is  the  atmospheric  density.  To 
determine  the  total  aerodynamic  torque,  the  individual  forces  must  be  integrated  as 

follows 


T aero  —  ^ R-cp-cg  ^  dF  (30) 

where  Rcp.<;g  is  the  vector  measured  from  the  center  of  mass  of  the  spacecraft  to  the  center 
of  the  individual  incremental  areas  (i.e.  the  distance  to  the  center  of  pressure  for  that 
individual  area). 
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The  aerodynamic  disturbance  torque  is  a  fixed  quantity  for  spacecraft 
which  are  oriented  toward  the  Earth,  but  varies  for  those  which  are  inertially  referenced 
such  as  this  spacecraft.  The  velocity  vector  for  the  spacecraft  can  be  expressed  in 
Earth-oriented  coordinates  as 


V,c  =  VXo 


(31) 


To  take  the  necessary  cross  products  for  determining  the  aerodynamic  torque,  this 
velocity  vector  must  be  expressed  in  body  coordinates  by  performing  the  two  required 
coordinate  transformations  using  the  Earth  to  sun  angles  and  the  sun  to  body  Euler 
angles.  Since  the  vector  r^p^g  for  each  individual  piece  of  the  satellite  is  known  in  body 
coordinates,  it  is  easy  to  simulate  the  cross  products  and  thus,  the  aerodynamic  torques 
for  each  individual  area  of  the  spacecraft.  These  individual  torques  are  then  added  to 
determine  the  overall  aerodynamic  torque. 

From  Reference  9,  the  worst  case  atmospheric  density  is  approximately 
2.0  X  10"'^  kg/m^.  Using  a  coefficient  of  drag  of  2.0,  orbital  velocity  of  7585.1  m/sec  and 
the  effective  surface  area  for  each  individual  satellite  panel  the  aerodynamic  force  on 
each  face  of  the  satellite  can  be  calculated  and  the  individual  cross  products  can  be 
determined.  The  individual,  time  varying  computations  are  performed  in  the  SIMULINK 
program  and  will  not  be  listed  here. 

b.  Solar  Pressure 

Radiation  incident  upon  the  body  of  the  spacecraft  produces  a  small 
disturbance  torque  about  the  center  of  mass  of  the  spacecraft.  The  surface  of  the  body  is 
subjected  to  a  radiation  pressure  (force  per  unit  area)  which  is  equal  to  the  vector 
difference  between  the  incident  and  reflected  solar  radiation  energy.  Since  the  intensity 
of  the  solar  radiation  is  a  function  of  distance  from  the  sxm,  spacecraft  altitude  with 
respect  to  the  Earth  does  not  affect  solar  radiation  pressure  for  Earth-orbiting  spacecraft. 
The  mean  momentum  flux  acting  on  a  surface  normal  is  determined  from 
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where  is  the  solar  constant  (1358  W/m^)  and  c  is  the  speed  of  light  (3.0  x  10*  m/sec). 
This  flux  value  is  constant  except  for  a  small  periodic  variation  due  to  the  eccentricity  of 
the  Earth's  orbit  about  the  sun  [Ref  6], 

For  most  applications,  the  solar  force  can  be  modeled  adequately  by 
assuming  the  incident  radiation  is  either  absorbed,  reflected  specularly,  reflected  diffusely, 
or  some  combination  of  the  above.  Specular  reflection  means  that  the  incident  vector  is 
reflected  at  the  same  angle  as  it  came  in  to  the  surface,  while  diffuse  reflection  means  that 
the  incident  vector  is  reflected  in  all  directions.  The  differential  force  imparted  by  the 
particles  which  are  totally  absorbed  is  given  by 


d¥abs  -  -PCaCOsQSdA  (33) 

where  is  the  absorption  coefficient,  S  is  the  unit  vector  from  the  spacecraft  to  the  sun 
and  0  is  the  angle  between  S  and  N,  the  unit  vector  normal  to  the  differential  surface.  This 
last  angle  must  be  less  than  90  degrees  for  a  pressure  force  to  exist  on  that  area.  The 
differential  force  for  the  specularly  reflected  radiation  is  given  by 

dPspec  =  PCs [S  -  2(cos  0N)  -  cos  0S]ci4  (34) 

where  C^is  the  coefficient  of  specular  reflection.  For  a  surface  which  reflects  diffusely,  all 
reflection  angles  must  be  considered  zmd  integrated  to  form  the  total  differential  force  for 
diffuse  reflection  as 


^/Frf,^=PO(-|cos0S-cos0S)ii4  (35) 

where  is  the  coefficient  of  diffuse  reflection.  If  all  types  of  reflection  play  a  part,  then 
the  total  differential  force  is  given  by 


total  =  -P  j[(l  -  ^5)8  +  2(Gcos0  +  }Q)N]coseci4 


(36) 


where  Cj+  C^+  Q  =  1  [Ref.  6], 

The  solar  torque  is  then  determined  by  taking  the  integral  of  all  cross 
products  between  each  differential  solar  force  and  its  individual  distance  offset  from  the 
center  of  mass  of  the  spacecraft  [Ref  6], 


Tcp-cg  '>^dF  total 


(37) 


Since  it  is  diflBcult  to  perform  this  integral,  it  is  not  unusual  to  break  down  the  spacecraft 
to  simple  geometric  elements  and  then  estimate  each  of  the  individual  torques  for  each 
element.  This  spacecraft  will  be  treated  in  a  similar  manner.  Since  only  the  solar  panel 
side  of  the  spacecraft  will  be  pointed  at  the  sun  during  normal  operations,  we  must 
consider  that  side  only  for  solar  pressure  (i.e.  the  angle  between  S  and  N  for  the  other 
sides  is  90  degrees  or  the  face  is  not  illuminated).  This  means  no  solar  torques  about  the 
Y  axis  will  be  generated.  The  surfaces  on  this  side  of  the  spacecraft  include  the  solar 
panels  and  the  sun-facing  side  of  any  monitoring  instruments  [Ref  6]. 

For  the  solar  arrays,  the  total  surface  area  is  1 .094  m^  and  it  is  assumed 
that  C,=  0.8,  C^=  0.1,  and  =  0.1.  The  center  of  pressure  for  the  array  assembly  is 
located  at  (x,  y,  z)  =  (0.0,  -9.5,  -14.9)  inches  or  (0.0,  -0.2413,  -0.3785)  meters  using  the 
body  coordinate  reference  listed  earlier.  Utilizing  Equations  (15)  and  (37),  this  will  cause 
a  torque  about  the  negative  X  and  Z  axes  equal  to 

Tspx  = -4. SlxlO-"^  N-m  (38) 

Tspz  =  -lAlxl0-^  N-m 

The  sun-monitoring  equipment  is  assumed  to  have  a  surface  area  of  0.061 
m^  and  it  is  assumed  that  C,=  0.05,  C,=  0.8,  and  =  0.15;  with  a  center  of  pressure 
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located  at  (x,  y,  z)  =  (0.0,  -9.5,  -32.25)  inches  or  (0.0,  -0.2413,  -0.8192)  meters.  This  will 
cause  a  torque  about  the  positive  X  axis  and  the  negative  Z  axis  equal  to 

T'eg'pwntt  —  1-92x10  '^N—m  (39) 

TggrpOTwte  =  ~1 -33x10  ^^N—m 

This  yields  a  total  solar  torque  for  each  axis  of 

Tsoiarx  =  -2.39x10“’^  N- m  (40) 

Tsolarz  -  -1-60x10"^  N-  m 

c.  Gravity  Gradient 

Any  non-symmetrical  body  in  orbit  will  experience  a  gravitational  torque 
because  of  the  variation  of  the  Earth's  gravitational  force  over  the  object.  This  is  due  to 
the  non-uniformity  of  the  inverse  square  gravitational  field.  To  develop  the  relationships 
for  the  gravity  gradient  torques  it  is  assumed  that  the  moments  of  inertia  are  known  and 
that  the  spacecraft  is  orbiting  a  non-realistic,  spherical  Earth.  Using  Reference  8,  the 
differential  force  acting  on  the  spacecraft  mass  element  dm  located  a  distance  R  from  the 
center  of  the  Earth  is  given  by 

=  (41) 

where  \ie  is  the  Earth's  gravitational  constant  and  R  is  the  distance  from  the  center  of  the 
Earth  to  the  individual  incremental  masses.  Using  Figure  10,  R  can  be  expressed  as 

R  =  Ro-r  (42) 

where  R^  is  the  distance  from  the  center  of  the  Earth  to  the  center  of  mass  of  the 
spacecraft  and  r  is  the  distance  between  the  spacecraft  center  of  mass  and  the  individual 
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incremental  masses.  Substituting  this  relationship  into  the  above  equation  and  using  the 
fact  that  r  /  Ro  «  1  to  expand  R‘^  in  a  Taylor  series  expansion  yields 


dF. 


^e(Ro-r)rfffi ..  >.r  Ro.. 

i\o 


(43) 


The  torque  about  the  center  of  mass  is  obtained  by  integrating  each  of  the 
individual  incremental  forces  and  noting  that 

J  r  Jw  =  0  (44) 

because  r  is  a  vector  measured  from  the  center  of  mass  of  the  spacecraft. 


Figure  10.  Gravity  gradient  spatial  relationships  [Ref.  8] 

The  final  result  is 

Tgg=  frxcflF  =  ^  f(rxRo)(r-Rc)c//«  (45) 

J  Ro  J 
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The  vector  Ro,  which  is  expressed  in  the  Earth-oriented  reference  frame  as  R,,  K,  can  be 
written  in  body  coordinates  by  performing  a  coordinate  transformation  from  the 
Earth-oriented  to  the  body  coordinate  system.  The  result  is 

AAA 

Ro  =  Ro{AX  l}Q(jy  +  IjoiJy  +  CWbod^  ('^6) 

where  the  expressions  for  A,  B  and  C  are  given  by 
A  =  5(1)65  [CGiSv]/  +  5(j)50Cv|/]  —  C(j)e5[50C<))] 

5  =  5(l)e5[CH/C(t)]  +  C(t)65[5'(l)]  (47) 

C  =  S^es[SQS^]l  -  5(l)C0Cv(/]  +  C<|)65[C0C(t)] 

where  C  is  the  cosine  of  the  angle  and  S  is  the  sine  of  the  angle.  Substituting  Equation 
(46)  into  Equation  (45)  and  integrating  yields  the  gravity  gradient  torques  about  each  axis 
in  body  coordinates  (ignoring  cross-product  terms)  as 


~  p3  Iyy)BC 

'^ggy  ~  ~  (48) 


’^gg^  ~  „3  -  lyy) AB 

JAo 


d.  Residual  Magnetic  Torque 

Through  the  use  of  non-ferrous  materials  and  proper  wiring 
configurations,  the  residual  magnetic  dipole  of  the  spacecraft  can  be  greatly  reduced  and, 
in  most  cases,  neglected.  There  can  be  no  specific  determination  of  the  residual  magnetic 
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dipole  until  the  spacecraft  is  actually  built  and  tested.  For  simulation  purposes,  residual 
magnetic  torques  will  be  ignored  in  this  application. 

e.  Total  Torque 

Initial  simulation  results  show  that  the  estimated  worst  case  torques  for 
each  axis  are  expected  to  be  approximately: 


roll 

2.0  X  10-' 

N-m 

pitch 

9.0  X  10-® 

N-m 

(49) 

yaw 

2.0  X  10-' 

N-m 

Since  the  yaw  and  roll  torques  are  expected  to  be  the  largest,  they  will  be  used  to  size  the 
ACS  equipment. 
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IV.  EQUATIONS  OF  MOTION  AND  ATTITUDE  CONTROL  LAWS 

A.  EQUATIONS  OF  MOTION 

From  Reference  8,  the  basic  moment  or  torque  equation  is 


Tc  =  ^crel  +  ®  X  He 


(50) 


where  is  given  for  a  system  with  a  momentum  storage  wheel  by  the  equation 


He  — /CD  +  h^y/Jee/  (51) 

where  I  is  the  inertia  matrix  of  the  satellite  and  co  contains  the  body  rotational  rates  of  the 
satellite  given  by 

C0;(  =  C0<j>  -  -  CO on(CQS\\i  +  5'([)iS'0C\j;) 

(Oy  =  Q-  5(1)1}/  -  CO  o„C(t)Ci|/  (52) 

CO2  =  50<i)  +  C(l)C0vi/  -  coo«(505v|;  -  5(()C0Cv|/) 


where  (j),  0  and  \|/  are  the  angles  used  in  the  3-1-2  Euler  angle  rotation  from  the 
sun-oriented  orbital  reference  frame  to  the  body  reference  frame.  The  term  (0o„  is  the 
angular  velocity  of  the  sun-oriented  reference  frame.  Assuming  the  time  rate  of  change 
of  Beta  to  be  zero,  it  is  given  by 


^on  —  0«  * 


cos^Gcoq 

cos^vsinp 


(53) 


Equation  (50)  can  be  written  in  components  as 
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Tx  —  (HyHz  ^ zHy 


Ty=Hy  +  (SizHx  -  (HxHz 


(54) 


Tz=Hz  +  (OxHy  -  (OyHx 

A  special  case  to  consider  occurs  when  the  rotating  axes  x,  y  and  z  coincide  with  the 
principle  inertia  axes  of  the  spacecraft.  Then  Equation  (51)  simplifies  to 


Hx  —  Ixx^x  +  hmfjggix 


Hy  —  lyydiy  4"  h 


^wheely 


(55) 


Hz  —  Izz^z  ^wheelz 


These  equations  can  be  substituted  back  into  Equation  (54),  noting  that  for  this 
application  there  is  only  wheel  momentum  about  the  y  axis,  to  form 


Tx  —  Ixx^x  +  ^y^zijzz  ~  lyy)  ~  (^zh^heely 


Ty  —  lyy&y  +  (Hx^zijxx  ~  Izz)  +  hyy/jggjy  (56) 

Tz  =  Izz^z  +  (^x^yifyy  ~  /xx)  +  ®x^»v/!ee/y 

These  equations  can  be  manipulated,  using  substitutions  from  above  and  noting  that  the 
wheel  momentum  is  oriented  along  the  negative  y  axis,  to  form  the  linearized  moment 
equations  for  each  major  axis 
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Tcx'^Tjx  ~  /xx($  ~  ®  onVj/  “  <0 0/1 +  (Jyy  ~  ^zz')(S^ o«^  +  ®  o/i<t*)  +  +  ®  o/i*|>) 

Tcy  +  Ti;ly  =  lyy^Q  —  G)  on^  ~  h  ^  (57) 

Tcz'^l'dz  ~  7zz(V  ®  on^  "I"  ®  0/1^)  ~  (7>y  ~  7xx)(®  o/i^  ~  ®  o/jV|/)  ~  /l((j)  ~  ©  o/iM^) 
where  T,.  constitutes  the  control  torques  and  constitutes  the  disturbance  torques. 

B.  CONTROL  FLOW 

Figures  1 1  and  12  present  the  control  flow  for  the  spacecraft  attitude  control 
system  and  Earth  sensor/Pitch  control. 


Figure  1 1 .  Spacecraft  attitude  control  flow 
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Figure  12.  Wheel/Pitch  control  flow 


C.  ATTITUDE  CONTROL  LAWS 

1.  Satellite  Capture  After  Booster  Separation 

Once  the  launch  vehicle  achieves  the  proper  Earth-oriented  configuration,  the 
satellite  will  be  separated.  Two  possible  separation  scenarios  will  be  considered  and 
simulated.  In  the  first  case,  the  momentum  wheel  will  be  spun  up  with  the  launch  vehicle 
providing  attitude  stabilization  until  the  wheel  reaches  optimum  speed.  Upon  separation, 
only  disturbance  torques  and  torques  generated  by  the  separation  mechanism  will  cause 
the  satellite  to  tumble.  For  the  second  case,  it  is  assumed  that  the  launch  vehicle  cannot 
provide  attitude  stabilization.  The  satellite  will  be  separated  from  the  launch  vehicle  and 
the  momentum  wheel  will  be  spxm  up  at  the  maximum  torque  rating.  This  procedure, 
along  with  any  torques  generated  by  the  separation  system  and  disturbance  torques,  will 
cause  the  satellite  to  spin  at  an  increasing  angular  rate  about  the  pitch  axis  until  the  wheel 
is  up  to  speed.  It  is  assumed  here  that  only  one  spring  failure  will  occur  among  the  four 
springs  in  the  separation  system.  This  single  spring  failure  will  cause  the  satellite  to 
tumble  if  the  wheel  is  not  up  to  speed,  and  will  cause  nutation  or  pitch  error  if  the  wheel 
is  up  to  speed.  All  tumbling  motion  needs  to  be  removed  so  the  satellite  can  be  stabilized 
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in  an  Earth-oriented  configuration  for  system  check  out  and  subsequent  transition  to  the 
sun-oriented  mode. 

Reference  2  suggests  the  use  of  the  following  magnetic  control  moments 

Mx  =  -KB:, 

My  =  -KBy  (58) 

Mz  =  -KBz 

where  K  is  the  control  gain  and  the  Bj  components  are  the  time  derivatives  of  the 
magnetic  field  relative  to  the  body  axes  of  the  satellite.  These  quantities  can  be  generated 
by  taking  the  derivative  of  the  magnetometer  outputs  [Ref  2]. 

This  theory  functions  on  the  premise  that  the  kinetic  energy  of  a  body,  due  to  the 
rotation  about  its  center  of  mass,  can  be  reduced;  which  eliminates  all  angular  motion 
between  the  satellite  and  the  magnetic  field.  This  results  in  a  terminal  motion  where  the 
satellite  rotates  twice  per  orbit  about  the  pitch  axis  (once  per  orbit  relative  to  the  Earth- 
oriented  orbital  reference  frame)  [Ref  2]. 

The  time  derivative  of  the  rotational  kinetic  energy  is  given  by 

^  =  T-o  (59) 

dt 

where  x  is  the  external  torque  acting  on  the  body  and  co  is  the  rotational  velocity  vector. 

A  torque  generated  by  the  magnetic  control  torque  rods  is  given  by 

x^ag  =  MxB  (60) 

Now  the  derivative  of  the  magnetic  field  can  be  estimated  by  (assuming  to»a><,) 
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B  ss  B  X  CO 


(61) 


Combining  the  above  equations,  yields 

^  =  B  •  M  (62) 

Implementing  the  magnetic  moment  equations  in  Equation  (58)  yields 

^  =  -KB'^  =  -A:|B  X  CO  P  (63) 

This  quantity  is  clearly  less  than  or  equal  to  zero  which  is  desired  to  remove  rotational 
energy.  Equation  (63)  shows  that  this  law  carmot  reduce  the  kinetic  energy  associated 
with  the  component  of  co  parallel  to  the  magnetic  field.  Fortunately,  the  orbital  motion  of 
the  satellite  will  ensure  that  the  direction  of  the  magnetic  field  relative  to  the  body  is 
continuously  changing  [Ref.  2]. 

Results  are  shown  in  Figures  13  through  15.  A  control  gain  of  1  x  10’  (N-m-sec^ 
/tesla^)  was  used  to  maintain  torque  rod  saturation  throughout  the  capture  procedure. 
Figure  13,  shows  the  results  of  a  simulation  where  the  wheel  has  already  been  spun  up 
and  the  initial  pointing  error  between  the  orbit  normal  and  the  spin  axis  is  20  degrees. 

This  was  done  to  show  that  the  system  can  remove  the  initial  nutation  due  to  the  inertial 
body  rotational  velocities  and  also  remove  the  pointing  error  between  the  orbit  normal 
and  spin  axis.  Visible  is  the  initial  chaotic  reaction  to  the  separation  where  the  rapidly 
changing  initial  response  is  due  to  nutation.  Once  all  body  rotational  rates  have  been 
neutralized,  the  system  begins  to  remove  the  pointing  error.  This  particulair  simulation 
took  approximately  2  orbits  to  remove  the  pointing  error  between  the  orbit  normal  and 
the  spin  axis. 
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orbit  normal/spin  axis  angular  error  in  degrees  vs  time 


Figure  13.  Initial  capture  pointing  error  to  orbit  normal.  Momentum  wheel 
spun  up  before  separation 

Figure  14,  displays  the  simulation  results  for  the  case  where  the  satellite 
momentum  wheel  has  not  been  spun  up  and  the  satellite  is  ejected  with  initial  rotational 
velocities  about  the  x  and  y  axes  of  approximately  five  degrees  per  second.  As  is  clearly 
apparent,  the  satellite  tumbles  through  several  revolutions  before  the  momentum  bias 
overcomes  the  body  angular  momentum  magnitude,  which  is  being  continuously 
decreased  via  the  control  law.  This  simulation  shows  that  the  satellite  went  through  30 
complete  rotations  before  it  settled  in  the  proper  configuration.  There  is,  however,  an 
inherent  problem  with  this  approach  to  capture  when  the  wheel  is  not  up  to  speed  at 
separation.  The  spin  axis  will  always  align  itself  with  the  orbit  normal,  but  there  is  the 
possibility  that  the  satellite  will  be  turned  about  the  yaw  axis  by  180  degrees  (momentum 
vector  pointed  in  opposite  direction  to  that  desired).  In  this  case,  ground  command  will 
be  necessary  to  process  the  momentum  vector  through  that  180  degree  swing  to  get  into 
the  proper  configuration.  Figure  15  is  an  expanded  view  of  Figure  14. 
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orbit  normal/spin  axis  angular  error  in  degrees  vs  time 


Figure  14.  Initial  capture  pointing  error  to  orbit  normal.  Momentum  wheel 
spun  up  after  separation. 


2.  Momentum  Wheel  Control 

Wheel  momentum  will  be  controlled  using  the  proportional/differential  control 


h  =  -KQ(TQQ  +  &)  (64) 

where  0  is  the  pitch  error  signal  provided  by  the  Earth  sensor  or  magnetometers  and 
Kq  and  xe  are  control  gains  [Ref.  8].  Substituting  Equation  (64)  into  Equation  (57)  yields 

Tcy  +  Tdy  —  lyy^ KqXqQ  +  KqQ  —  Iyy(bon  (65) 

Ignoring  the  wheel  desaturation  control  input  T<,  for  now  and  considering  the  last  term  on 
the  right  hand  side  to  be  part  of  the  disturbance  torques  (this  term  is  on  the  order  of  10'^  in 
magnitude  or  it  can  be  canceled  by  feeding  forward  an  estimate  of  the  signal)  yields 
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orbit  normal/spin  axis  angular  error  In  degrees  vs  time 


Figure  15.  Exploded  view  of  Figure  14  showing  final  pointing  angle. 

Tdynew  —  +  X^eT00  +  KqQ  (66) 

Transforming  to  the  Laplace  domain 
Qt'y)  _ _ 

r dynewi^S)  J^O(5)+^0T0//^J0(j)+/ro//^0(5) 

This  is  the  transfer  fimction  for  a  second  order  damped  oscillator  with  natural  frequency 


and  damping  ratio 


The  design  point  is  for  a  critically  damped  system  where  ^=1  and  then 
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(70) 

If  Kg  is  selected  as  0.36  kg-m^/sec^  and  Xg  is  selected  as  10  seconds,  the  natural  frequency 
is  0.19  radians/second  with  a  settling  time  of  approximately  40  seconds,  assuming 
settling  time  is  defined  as  four  times  the  time  constant  Xg. 

System  response  to  an  initial  pitch  error  of  five  degrees  is  shown  in  Figure  16. 

The  system  is  reacting  as  a  critically  damped  oscillator  with  a  settling  time  of 
approximately  40  seconds,  as  designed. 


theta  angle  response  to  initial  pointing  error 


Figure  16.  Momentum  wheel  control  response  to  initial  pitch  error. 


An  interesting  pitch  response  occurs  when  Earth  sensor  data  is  unavailable  and 
only  magnetometer  pitch  error  is  fed  back  for  control.  Since  pitch  rate  is  no  longer 
available,  Equation  (67)  becomes  an  undamped  oscillator  which  will  grow  over  time  vnth 
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pitch  disturbance  torque  input.  Fortunately,  Earth  sensor  data  is  unavailable  for  only  a 
short  period  of  time,  so  the  pitch  response  does  not  grow  too  large  as  in  Figure  17.  If  this 
oscillation  becomes  a  problem  for  the  overall  system  response,  then  some  rough 
estimation  of  pitch  rate  will  need  to  be  generated  to  provide  some  damping. 


pitch  response  after  magnetometer  data  feedback 


50  55  60  65  70  75  80 

time  in  minutes 


Figure  17.  Pitch  response  with  magnetometer  pitch  angle  feedback  only 


3.  Momentum  Wheel  Desaturation  Control 

If  a  secular  disturbance  torque  is  experienced  about  the  pitch  axis,  momentum 
wheel  speed  will  grow  without  bound  to  counteract  the  disturbance  torque  [Ref.  8].  Once 
the  momentum  wheel's  speed  has  deviated  far  enough  away  from  its  nominal  speed  (3000 
rpm),  something  must  be  done  to  return  the  wheel  to  the  correct  speed.  That  is  the 
purpose  of  the  momentum  wheel  desaturation  control.  The  portions  of  the  magnetic 
control  moments  for  wheel  desaturation  are  given  by 


47 


Mx  =  -KuBz^h 


My  =  0 


(71) 


Mz  =  KuBxAh 


where  Ah  is  the  overall  change  in  the  wheel's  momentum  from  its  nominal  value  and  K^is 
the  control  gain  for  desaturation.  After  the  cross  product  is  taken  with  the  magnetic  field, 
the  control  torques  are  calculated  to  be 

Tx  =  —KuBxByAh 

Ty  =  Ku(.Bl+BhAh  (72) 

Tz  =  -KuByBzAh 

Of  note,  this  control  law  generates  the  required  torque  about  the  y  axis  for  wheel 
desaturation,  but  it  also  creates  disturbances  about  each  of  the  other  axes  which  must  be 
accounted  for.  Also,  the  control  gain,  ,  must  be  a  negative  number.  This  can  be  seen 
by  the  need  for  a  negative  y  axis  torque  if  Ah  is  a  positive  quantity  and  vice  versa.  The 
control  gain  K„is  sized  to  provide  torque  rod  saturation  when  Ah  is  ten  percent  of  nominal 
wheel  momentum.  This  leads  to  a  control  gain  of  -1.538  x  10® (sec  '  -  tesla'^). 

Momentum  wheel  desaturation  control  will  be  shut  off  when  sun  sensor  data  is 
unavailable  because  the  control  law  creates  disturbances  for  the  x  and  z  axes  which  cannot 
be  dissipated  while  there  is  no  sun  sensor  data. 

4.  Control  for  Transition  from  Earth-oriented  to  Sun-oriented  Attitude 

The  control  law  for  sun-oriented  control  was  developed  first  with  the  hope  that  it 
could  be  used  without  modification  to  accomplish  the  transition  from  Earth-oriented/ 
checkout  operations  to  sun-oriented  operations.  The  orbit  configuration  will  be  adjusted 
in  order  to  vary  the  Earth  to  sun  angles  (t)„  and  and  test  the  control  law  at  several 
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values  of  Beta  angle.  It  was  not  known  whether  the  sun-oriented  control  law  could 
handle  such  large  angular  errors  since  the  control  law  was  based  on  small  angle 
assumptions). 

Reorientation  simulations  were  only  performed  for  a  1 :30  orbit  and  a  3:00  orbit. 

It  was  assumed  that  if  the  reorientation  worked  for  these  initial  conditions,  then  it  would 
work  for  the  4:30  and  6:00  orbit  because  the  transition  angles  will  be  smaller  and 
disturbance  torques  will  be  smaller  due  to  the  lower  Earth  to  sim  angular  values. 

Figure  18  displays  the  results  of  a  reorientation  simulation  from  an  initial  sun 
pointing  error  of  approximately  72  degrees.  This  constitutes  a  reorientation  for  a  1 :30 
am/pm  orbit  with  average  Beta  angle  of  18  degrees.  It  takes  several  orbits  to  perform  the 
reorientation  as  expected.  Figure  19  shows  a  reorientation  simulation  for  a  3:00  am/pm 
orbit.  The  raised  dotted  lines  at  the  bottom  of  the  plots  represent  the  time  periods  when 
the  satellite  is  not  in  eclipse.  If  the  value  of  the  dotted  line  is  one,  then  roll/yaw  attitude 
sensing  data  is  available.  The  eclipse  control  also  serves  to  cut  off  momentum  wheel 
dumping  control  when  the  satellite  is  in  eclipse. 

5.  Sun-oriented  Attitude  Control 

The  pitch  axis  control  torque  will  be  generated  as  discussed  in  the  momentum 
wheel  and  wheel  desaturation  control  sections.  The  magnetic  moment  equations  are 
given  as 


Mx  =  -KuBzAh 


My  =  K\Bx^  +  ~  Bx^)  +  K^Bz^f 


(73) 


Mz  =  KuBxMi 


where  K,  and  K3  are  used  for  precession  control,  while  K2  is  used  to  control  nutation. 
Reference  1 0  used  the  control  law 
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K4y  —  K\ Bx^  +  K2By 


(74) 


For  small  angles,  since  no  measurement  for  yaw  existed;  By  can  be  approximated  by 


By  =  -  C0o(^  +  2(1)5;,) 


(75) 


In  the  orbital  reference  frame,  By  is  constant  for  an  Earth-pointing  satellite.  For  a 
sun-pointing  satellite,  it  is  not  constant,  so  the  control  law  from  Reference  10  should  not 
be  used.  Instead,  the  y  axis  torque  rod  control  moment  is  selected  as  the  second  equation 
of  Equation  (73),  which  is  quite  similar  to  Reference  10.  The  additional  term  has 

been  added  because  there  is  a  yaw  error  measurement  which  is  not  available  when  using  a 
horizon  sensor,  which  is  the  standard  sensor  for  Earth-pointing  satellites. 


Figure  18.  Earth  to  sun  reorientation  (1:30  orbit) 
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When  the  M  X  B  cross  product  is  performed  the  following  control  torques  are  generated 
Tx  =  K\BxBz^  +  K2{B'z^  ~  BxBz'^')  ~  KuBxBy^h 

Ty^Ku{Bl+Bl)Lh  (76) 

Tz  —  —K\Bx^  ~ K2{BxBz^  ~~ -fixy)  —  K.'^BxBz'^ ~ KuByBz^h 

When  these  control  torques  are  substituted  into  Equation  (57)  above,  the  total  torque 
equations  become 

Tdx—^xx^~  K.2B^^  —  [(/xx  —Iyy  +  /zz)®  on  +  h—  K2BxBz]'^  +  KuBxBz^ 

+{(Jyy  -  7zz)®  \n+h(£ton-K\ BxBz]^  -  [/xxtO  on+K^Bz]'^ 

Tdy  —  lyy^  +/l9T60  +  KqQ  —  KuiBx  +  Bz)^  ~  lyy^ on  (77) 
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Tdz  ~  /zzV  [(-^zz  ^yy  ~  h  +  IC2BxBz^^  ~  K2Bx'^  +  KuByBzAh 

"^[Izz^ on  "I"  K\Bx^^  +  [(/xv  ~  Ixx)^on  +  h(£)on  +  K.^BxBz^'^ 


Figures  20  and  21  show  the  reason  for  selecting  the  control  laws  as  listed.  In  Figure  20,  a 
positive  (j)  error  is  shown  which  requires  a  change  in  momentum  in  the  positive  z 
direction.  This  can  be  generated  by  the  (j))  portion  of  the  z  axis  magnetic  torque  in 

Equation  (76).  If  K,  is  selected  as  a  negative  number,  this  generates  a  positive  torque  as 
required  to  achieve  the  desired  change  in  momentum.  Figure  21  shows  what  happens 
when  there  is  a  negative  \\f  error.  This  error  requires  a  change  in  momentum  in  the 
positive  X  direction  which  is  generated  by  the  (KjB/  H/)  portion  of  the  x  axis  control  law 
in  Equation  (76).  Of  note,  the  portions  of  the  x  and  z  axes  control  laws  which  contain 
B^B^  create  disturbance  torques  about  their  respective  axis.  However,  over  the  course  of 
one  orbit,  these  disturbances  should  average  out  to  zero.  This  averaging  could  be 
precluded  by  the  effects  of  eclipse  which  would  cut  off  attitude  data  for  a  portion  of  the 
orbit. 


momentum  vector  (h)  into  page 

X 

A  h  located  in  the  plane 
of  the  page 

/  ^ 

1  _ 

7 

Figure  20.  Generation  of  required  control  torques;  z  axis 
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momentum  vector  (h)  into  page 


X 


The  roll  and  yaw  equations  can  be  expressed  in  matrix  form  as 

T(ix~KuBxBzAh  ^  A  C  (78) 

Tdz-KuByBzAh  J  "  [  5  £>  J[  miCj)  J  ^  ^ 

where 

A  =  IxxS^  ~  K2B'zS  +  \(Jyy  ~  Izz)^  on  +  hOi  on  ~  K\BxBz'\ 

B  =  [(Jzz  ~  lyy  +  /xx)©  on~h  +  K2BxBz'\s  +  \Izz^  o«  +  -^1  Bx  ]  ("79) 

C  =  ~[(/xx  “  lyy  +  fzz)®  on'^h  —  K2BxBz\s  —  [/xx®  on  +  K^Bz ] 

D  =  IzzS^  ~  K2BxS  +  \jjyy  ~  fxx)®  o«  +  h(£>on  +  K^BxBz'l 

Using  methods  introduced  in  Reference  10,  the  nutation  and  orbital  time  constants  of  the 
system  can  be  estimated  by 

^ _ 4hI„Izz _ 

Bo[(4/xt+/K)/l^2~^xr^ir(4/k3+/L  ] )] 
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(80) 


4/70)0 

T^Ol  =  ““1 - 

nBl{Kx+\OK^) 

4/70)0 

Xo2  = - ^ - 

where  Bj,  is  M/R^  of  Equation  (26)  and  the  time  constants  are  expressed  in  orbits.  The 
values  of  K,  and  K3  were  selected  to  saturate  the  torque  rods  for  an  angular  error  of  one 
degree.  This  leads  to  a  selection  of  the  control  gains  as 

K,  =  -  9.555  X  10'  Amp-mVTesla 

K2  =  -  4.777  X  10*  Amp-m'  -sec/Tesla  (81) 

K3  =  -  2.389  X  10'  Amp-m'/Tesla 

which  results  in  time  constant  values  of 
Tn  =  13.711  seconds 

T,  =  6.848  minutes  (82) 

Tj  =  5.326  minutes 

These  time  constants  only  apply  under  conditions  of  no  torque  rod  saturation  and 
continuous  control  (i.e.  no  eclipse  shut  off  of  attitude  sensing  data).  Since  the  gains  are 
sized  to  create  torque  rod  saturation  above  one  degree  pointing  error,  the  simulated  time 
constants  will  be  slightly  longer.  Of  note,  K,  AND  K3  are  destabilizing  for  nutation. 

This  can  be  managed  by  ensuring  that  the  following  relationship  is  satisfied 

(4K3+K,)<[(4I„+IJhKJ/I„I„  (83) 

Figures  22  and  23  show  simulation  results  which  can  be  used  to  estimate  system 
time  constants.  Figure  22  shows  nutation  and  Figure  23  shows  the  orbital  time  constants 
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Using  these  figures  to  estimate  the  time  constants  yields  a  nutation  time  constant  of 
approximately  21  seconds,  and  an  orbital  time  constant  of  approximately  6.5  minutes. 
These  estimations  are  close  to  the  calculated  values. 


pointing  error  in  degrees  vs  time;  sun-oriented  control 


Figure  22.  Nutation  time  constant  response 

Figures  24  through  34  show  simulation  results  for  sun-oriented  control  in  different 
orbit  configurations  (1 :30,  3:00, 4:30  and  6:00  orbits).  Figures  24  through  27  do  not 
include  sensor  noise  inputs  or  cross  product  of  inertia  effects.  Figures  28  through  3 1  do 
include  these  effects  with  a  sensor  noise  gain  of  0.07  degrees  which  is  much  higher  than 
those  typically  found  in  commercial  sensors.  Figure  28  and  Figures  32  through  34 
include  results  for  a  1 :30  orbit  with  noise  and  cross  product  of  inertias  at  different  Earth 
positions  around  the  sun.  Figure  28  simulates  Summer  Solstice,  Figure  32  simulates 
Vernal  Equinox,  Figure  33  simulates  Autumnal  Equinox,  and  Figure  34  simulates  Winter 
Solstice.  All  eleven  simulations  use  the  following  initial  conditions 


(j)o  =  1  degree 

co,o  =  0.00075 

radians  /  second 

00  =  2  degree 

®y0  ®on0 

radians  /  second  (see  Equation  (53)) 

vj/o  =  -0.5  degree 

©30  =  -0.00003 

radians  /  second 
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sun-oriented  control;  orbital  time  constants 


Figure  23.  Orbital  time  constant  response 


Figiue  24.  Sun  oriented  control;  1 :30  orbit;  no  noise  input 
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pointing  error  in  degrees  vs  time;  sun-oriented  control 
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Figure  28.  Sun  oriented  control;  1:30  orbit;  with  noise  input 
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pointing  error  in  degrees  vs  time;  sun-oriented  control 


Figure  31.  Sun  oriented  control;  6:00  orbit;  with  noise  input 


pointing  error  in  degrees  vs  time;  sun-oriented  control 


Figure  32.  Sun  oriented  control;  with  noise  input.  Vernal  Equinox 
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pointing  error  in  degrees  vs  time;  sun-oriented  control 


Figure  33.  Sun  oriented  control;  with  noise  input.  Autumnal  Equinox 


Figure  34.  Sun  oriented  control;  with  noise  input.  Winter  Solstice 
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V.  CONCLUSIONS  AND  RECOMMENDATIONS 


A.  CONCLUSIONS 

Equations  of  motion  for  a  sun-oriented  inertial  reference  frame  were  developed. 
Transformation  to  a  body-oriented  reference  frame  was  accomplished  using  a  3-1-2  Euler 
axis  rotation  sequence  with  small  angle  approximations. 

Magnetic  attitude  control  laws  for  roll  and  yaw  control  were  developed  by 
modifying  control  laws  which  have  been  used  in  the  past  for  Earth-oriented  operation. 
Pitch  control  was  accomplished  by  varying  momentum  wheel  speed  using  a  proportional  / 
derivative  feedback  controller.  The  Earth-orientation  control  laws  work  well  for  this 
sun-oriented  control  application. 

Simulations  were  performed  to  test  capture  control  of  the  satellite  once  it  is 
separated  from  the  launch  vehicle.  The  satellite  was  separated  wtih  and  without  the 
momentum  wheel  being  spun  up  to  nominal  operating  speed  to  test  system  performance. 
The  control  law  is  capable  of  removing  all  rotational  kinetic  energy  and  aligning  the 
momentum  vector  along  the  orbit  normal.  The  possibility  exists  for  the  momentum  vector 
to  be  aligned  away  from  the  sun  (180  degree  misdirection)  when  the  satellite  is  separated 
before  the  momentum  wheel  is  spun  up. 

Reorientation  from  Earth-oriented  /  system  check-out  operations  to  a  sun-oriented 
mode  was  performed  using  the  same  magnetic  attitude  control  law  as  used  for  the 
sun-oriented  roll  and  yaw  control.  The  control  law  does  not  break  down  for  low  Beta 
angle  orbit  applications. 

Sun-oriented  control  maintains  a  pointing  accuracy  better  than  0.2  degrees  for  the 
worst  case,  1:30  orbit  configuration.  Given  this,  momentum  wheel  size  could  probably  be 
cut  in  half  to  maintain  the  desired  pointing  accuracy  within  0.5  degrees.  This  would 
reduce  weight  and  cost  while  also  reducing  reorientation  times  from  the  Earth-oriented  to 
the  sun-oriented  mode.  A  second  Earth  sensor  is  not  required  to  maintain  pointing 
accuracy  during  eclipse  in  the  current  configuration. 
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Noise  and  cross  product  of  inertia  inputs  do  not  have  a  noticeable  affect  on 
pointing  accuracy  output. 

B.  RECOMMENDATIONS 

Future  study  could  be  done  which  implements  a  Kalman  Filter  for  noise 
attenuation  in  the  sensors.  The  filter  could  also  be  used  to  estimate  pitch  rate  from 
magnetometer  data  when  Earth  sensor  data  is  unavailable.  This  would  avoid  the 
undamped  oscillations  experienced  without  pitch  rate  feedback. 

Study  could  be  performed  using  a  6:00  orbit  (continuous  sun  sensor  input)  without 
a  momentum  bias  to  test  the  ability  of  the  torque  rods  to  control  the  spacecraft  by 
themselves.  Torque  rod  size  could  also  be  modified  to  see  how  small  the  torque  rods  can 
be  before  they  are  unable  to  control  the  satellite. 

The  current  SIMULINK  program  could  be  modified  in  order  to  speed  up 
simulation  times  by  removing  several  of  the  function  blocks  and  constructing  "hard  wired" 
representations. 
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APPEIVDIX:  SIMULATION  BLOCKS  AND  MATLAB  CODE 

MATLAB/SIMULINK  block  diagrams  are  included  here  for  the  program  used  to 
simulate  all  satellite  control  sequences.  Computer  code  for  MATLAB  control  programs 
are  also  listed.  The  main  SIMULINK  program  block  is  shown  in  Figure  A1 .  The  beta 
angle  and  Earth  to  Sun  angle  gemorater  block  creates  beta  angle  and  Earth  to  sun  angles 
in  accordance  with  Chapter  II  equations.  The  B-field  block  generates  the  Earth's 
magnetic  dipole  field  in  Earth-oriented  orbital  coordinates.  This  field  is  then  transformed 
to  body  coordinates  by  performing  the  two  coordinate  transformations;  3 -axis  orbit  to  sun 
frame  and  3 -axis  sun  to  body  frame.  The  magnetometer  block  acts  as  the  magnetic  field 
sensor  and  time  derivative  estimator  for  the  rate  of  change  of  the  magnetic  field  with 
respect  to  the  spacecraft  body.  These  outputs  are  then  fed  to  the  Magnetic  torque  rod 
control  along  with  change  in  wheel  momentum  (Ah),  spacecraft  attitude  and  spacecraft 
attitude  rate.  Magnetic  moments  are  then  created  and  fed  to  the  torque  rod  block  where 
saturation  limiting  is  performed.  The  saturation  limited  magnetic  moments  are  then  used 
by  the  cross  product  block  M  x  B  to  calculate  magnetic  control  torque  which  is  the  cross 
product  between  the  magnetic  moments  and  the  magnetic  field  expressed  in  body 
coordinates.  These  torques  along  with  the  disturbance  torques  are  sent  to  a  torque 
summer  which  feeds  the  overall  external  torque  input  to  the  spacecraft  dynamics  block 
where  the  body  rotational  rates  (body  rates)  are  calculated.  The  momentum  from  the 
momentum  wheel  is  also  input  here  to  satisfy  the  equations  of  Chapter  IV.  The  body 
rates  are  then  input  to  the  kinematics  section  to  calculate  the  Euler  angle  rates  which  are 
integrated  to  form  the  Euler  angles  (j),  0  and  \\i.  The  attitude  sensor  block  passes  through 
attitude  information  with  sensor  noise  input  when  sensor  data  is  available.  At  all  other 
times,  is  passes  an  error  signal  of  zero  to  simulate  the  loss  of  attitude  sensing  information. 
Individual  output  blocks  for  each  variable  to  be  plotted  are  included  to  pass  data  to  the 
MALAB  workspace  for  subsequent  processing.  Each  individual  block  is  displayed  on  the 

following  pages. 
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Figure  A1 .  Overall  SIMULINK  block  diagram 
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Figure  A2.  Earth's  magnetic  field  in  Earth-oriented  coordinates 


Figure  A3.  Beta  and  Earth  to  sun  angle  generation 
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Figure  A4.  Beta  angle  generation 


Figure  A5.  Coordinate  transformation  from  Earth  orientation  to  Sun  orientation 
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Figure  A6.  Coordinate  transformation  from  Sun  orientation  to  body  orientation 
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Figure  A8.  Magnetic  field  derivative  approximation  by  cross  product  with  Euler  rates 


Figure  A9.  Sun-oriented  magnetic  torque  rod  control 
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Figure  A14.  Attitude  sensors.  Earth  and  Sun  sensors  with  noise  input 
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Figure  A15.  Momentum  wheel  control 


Figure  A 16.  Momentum  wheel  saturation  controls 
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F igure  A 1 7 .  Spacecraft  Dynamics 


Figure  A1 8.  Spacecraft  Kinematics 
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Figure  A 19.  Disturbance  torque  generation 


Figure  A20.  Gravity  gradient  torque  generator 
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Figure  All.  Aerodynamic  disturbance  torque  generator 


Figure  All.  Generation  of  satellite  velocity  unit  vector  in  body  coordinates 
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Figure  A23.  Computation  for  Aero  disturbance  torque  for  each  area 
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MATLAB  CODE 


%this  program  is  the  controller  for  a  simulink  file  projconS  which  will  simulate 
%the  equations  of  the  attitude  determination  and  control  system 

%The  program  will  set  up  several  variables  and  initial  conditions  and  will 
%also  control  data  display 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

printcon=0; 

control=l ;  %set  to  1  to  call  simulink  program 

mode=0; 

lawcontrol=l ; 

if  control=l, 

clear 

control=l; 
lawcontrol=l ; 
mode=0; 
printcon=0; 
end 

orient  tall 

%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%variable,  constants  and  initial  conditions 

maxrows=[80000,8];  %sets  max  #  of  data  points 
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stepsizemin=.09; 

stepsizeinax=.09; 

starttime=0; 

stoptime=10*5760; 

eclipsecon=0; 

rearth=6378.135; 

alt=550; 

earthrotvel=(3 60/24)*  ( 1/3 600)*  (pi/ 180); 

wo=sqrt(3 9860 1 .2/(rearth+alt)^3); 

rad/sec 

omega0=wo; 

inc=97.59306*pi/180; 

nuO=wo*starttime*pi/l  80; 

node 

epsilon=23.442*pi/l  80; 

% 

%nioments  of  inertia 

ixx=8.363; 

iyy=9.284; 

izz=4.012; 

ixy=-.lll; 


%set  step  size  for  simulation 


%1  means  no  eclipse 
%radius  of  Earth  in  km 
%satellite  altitude  in  km 
%spin  rate  of  earth  in  rad/sec 
%orbital  rotational  velocity  in 

%orbit  inclination  in  degrees 
%angle  around  orbit  from  ascending 

%earth  tilt  wrt  ecliptic 

%kg-m'^2 
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ixz=.013; 


iyz=.048; 

ii=[ixx  -ixy  -ixz;-ixy  iyy  -iyz;-ixz  -iyz  izz]; 

invl=inv(ii);  %inverse  of  inertia  matrix  for  dyn 

% 

%beta  angle  generation 

uoffset=88;  %set  inital  Earth  position  around  sun 

betaw0=(uoffset+22.5)*pi/l  80; 

betawdot=-9.9639*cos(inc)*(l/24)*(l/3600)*(pi/180)*(rearth/(rearth+alt))^3.5; 
betau0=(uoffset+0)*pi/l  80; 
betaudot=(360/365.242)*(l/24)*(l/3600)*(pi/180); 

betaO=asin(cos(inc)*sin(betauO)*sin(epsilon)+sin(inc)*sin(betawO)*cos(betauO)... 
-cos(epsilon)  *  cos(betawO)  *  sin(betauO)  *  sin(inc)) ; 

% 

%solve  for  nusim  for  earth  to  sun  angle  computations 

A=acos(-cos(pi-inc)  *  cos(epsilon)+sin(pi-inc)*  sin(epsilon)  *  cos(betawO)) ; 

phro=asin(sin(pi-inc)*sin(betaw0)/(sin(A)+.0000000001))-betau0; 

B=asin(sin(epsilon)*sin(betaw0)/(sin(A)+.00000000001)); 

bbb=tan(betaO); 

ccc=tan((pi/2)-A); 

eee=asin(bbb  *  ccc) ; 
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nusun=B-(eee); 

nuesO=nuO+nusun; 

% 

%solve  for  theta  earth  to  sun  at  time  0 
thetaesO=atan2(tan(nuesO),sin(betaO)); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%disturbance  torques 

%center  of  mass  listings  in  meters 

cgx=  01*.0254; 

cgy=-.23*.0254; 

cg2r=-l  7.84’'‘.0254; 

% 

%disturbance  torque  control 

%aero  cutoffs 

kal=l; 

ka2=l; 

ka3=l; 

ka4=l;  . 

ka5=l ; 

ka6=l; 

ka7=l; 
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%axis  cutoff  by  type  of  disturbance  torque 


ggxcon=l; 

soIarxcon=l; 
aeroxcon=l; 
ggycon=l; 
solary  con=l; 
aeroycon=l; 
ggzcon=l; 
solarzcon=l ; 
aerozcon=l; 

% 

%axis  torque  cutoff 
distxcon=l; 
distycon=l; 
distzcon=l; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
%sensor  noise 

earthsenscon=0*(.07*pi/l  80);  %set  nonzero  for  noise  input 

rollsunsenscon=0*(.07*pi/l  80); 
yawsunsenscon=0*(.07*pi/l  80); 

% 
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%b-field  generation  initial  conditions 

alfaO=nuO; 

mu0=0; 

% 

%wheel  control  gains  and  momentum  bias 
taup=10; 

Kp=3.6e-1; 

wheelbias=9.75;  %wheel  inertia  is  .0310352  kg-m^2 

wheelh0=0; 

% 

%torque  rod  control 
Kl=(-5*1.911e7); 

K2=-5*9.554e7; 

K3=(K1)/4; 

Ku=-l*(1.538e6); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%kinematics  initial  conditions 

phi0=(l*pi/180); 

theta0=(2*pi/180); 

ksi0=(-.5*pi/180); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
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%spacecraft  dynamics  initial  conditions 
wx0=0; 

wy0=((cos(thetaes0)^2)*(-wo)/((cos(nues0)^2)*sin(beta0))); 

wz0=0; 

% 

HxO=ixx*wxO; 

HyO=iyy*wy  0-9.75 ; 

Hz0=izz*wz0; 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%i 
if  control=l, 

% 

%call  simulink  program 

fmalprS 

pause 

end 

% 

wheelspeed=wheelmomentum*60/(2*pi*  .03 1 03  52); 
%percentspeed=100*abs(abs(wheelspeed)-3000)/300; 

%percentwheeltorque=  1 00  *  abs(hdotwheel)/.02 5 ; 

%percentdeltamomentum=  1 00*  abs(wheeldeltamomentum/. 9749996); 

% 
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%phisqrd=phi.  *phi; 

%thetasqrd=theta.*theta; 

%ksisqrd=ksi.*ksi; 

pointerror=sqrt((phi.*phi)+(ksi.*ksi)); 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

%plot  results 

clg 

if  control=101, 

subplot(3 1  l),plot(time/60,bxorb) 
title('bfild  in  Earth-oriented  frame') 
xlabel('time  in  minutes') 
ylabel('bfield  strength  x  in  Earth-orbit  frame') 
grid 

subplot(3 1 2),plot(time/ 60, byorb) 
xlabel('time  in  minutes') 
ylabel('bfield  strength  y  in  Earth-orbit  frame') 
grid 

subplot(3 1 3),plot(time/60,bzorb) 
xlabel('time  in  minutes') 
ylabel('bfield  strength  z  in  Earth-orbit  frame') 
grid 
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pause 


ifprintcon=l, 

print 

end 

% 

subplot(3 1 1  ),plot(time/60,bxsun) 
title('bfield  in  Sun-oriented  frame') 
xlabel('time  in  minutes') 
ylabel('bfield  strength  x  in  sun-oriented  frame') 
grid 

subplot(3 1 2),plot(time/60,bysun) 
xlabel('time  in  minutes') 
ylabel('bfield  strength  y  in  sun-oriented  frame') 
grid 

subplot(3 1 3),plot(time/60,bzsun) 

xlabel('time  in  minutes') 

ylabel('bfield  strength  z  in  sun-oriented  frame') 

grid 

pause 

if  printcon=l, 
print 
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end 


% 

subplot(3 1  l),plot(time/60,bxbody) 
title('bfield  in  body  frame') 
xlabel('time  in  minutes') 
ylabel('bfield  strength  x  in  body  frame') 
grid 

subplot(3 1 2),plot(time/60, bybody) 
xlabel('time  in  minutes') 
ylabel('bfield  strength  y  in  body  frame') 
grid 

subplot(3 1 3),plot(time/60,bzbody) 
xlabel('time  in  minutes') 
ylabel('bfield  strength  z  in  body  frame') 
grid 
pause 

if  printcon— 1, 

print 

end 

end 
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%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

subpIot(3 1  l),plot(time/60,o2sphi*  1 80/pi) 

title('orb  to  sun  angles  vs  time') 

xlabel('time  in  minutes') 

ylabel('o2s  phi') 

grid 

subplot(3 1 2),plot(time/60,o2stheta*  1 80/pi) 
xlabel('time  in  minutes') 
ylabel('o2s  theta') 
grid 

subplot(3 1 3),plot(time/60,o2sksi*  1 80/pi) 

xlabel('time  in  minutes') 

ylabel('o2s  ksi') 

grid 

pause 

if  printcon=l, 

print 

end 

% 


clg 

subplot(3 1 1  ),plot(time/60,o2sthetadot) 


xlabel('time  in  minutes') 
ylabel('o2s  thetadot') 
grid 

subplot(3 1 2),plot(time/60,betaangle*  1 80/pi) 
xlabel('time  in  minutes') 
ylabel('betaangle') 
grid 

subplot(3 1 3),plot(time/60,woddotfunet) 
xlabel('time  in  minutes') 
ylabel('omegaOnew  dot  rad/see''^2') 
grid 
pause 

if  printcon=l, 
print 
end 
clg 
% 

if  control=101 

subplot(3 1  l),plot(time/60,velunitvect(:,l)) 
title('velocity  unit  vector  components  in  body  frame') 
xlabel('time  in  minutes') 
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ylabel('velocity  in  x  body') 
grid 

subplot(3 1 2),plot(time/60,velunitvect(:,2)) 
xlabel('time  in  minutes') 
ylabel('velocity  in  y  body') 
grid 

subpIot(3 1 3),plot(time/60,velunitvect(:,3)) 

xlabel('time  in  minutes') 

yIabel('velocity  in  z  body') 

grid 

pause 

if  printcon=l, 

print 

end 

end 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 

subplot(3 1 1  ),plot(time/60,phi*  1 80/pi) 

title('euler  angle  in  degrees') 

xlabel('time  in  minutes') 

ylabel('roll  angle  in  degrees') 

grid 
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subplot(3 1 2),plot(time/60, theta*  1 80/pi) 
xlabel('time  in  minutes') 
ylabel('pitch  angle  in  degrees') 
grid 

subplot(3 1 3),plot(time/60,ksi*  1 80/pi) 
xlabel('time  in  minutes') 
ylabel('yaw  angle  in  degrees') 
grid 
pause 

if  printcon=l, 
print 
end 
% 

subplot(3 1 1  ),plot(time/60,phidot*  1 80/pi) 
title('euler  angle  rates  in  degrees  per  second') 
xlabel('time  in  minutes') 
ylabel('roll  angle  rate  in  deg/sec') 
grid 

subplot(3 12),plot(time/60,thetadot*  1 80/pi) 
xlabel('time  in  minutes') 
ylabelCpitch  angle  rate  in  degrees') 
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grid 


subplot(3 1 3),plot(time/60,ksidot*  1 80/pi) 

xlabel('time  in  minutes') 

ylabel('yaw  angle  rate  in  degrees') 

grid 

pause 

if  printcon=l, 
print 
end 
% 

if  control=101, 

subplot(3 1  l),plot(time/60,ggx) 

titleCgrav  grad  torques  vs  time') 

xlabel('time  in  minutes') 

ylabel('grav  grad  x') 

grid 

subplot(3 1 2),plot(time/60,ggy) 
xlabel('time  in  minutes') 
ylabel('grav  grad  y') 
grid 

subplot(3 1 3),plot(time/60,ggz) 
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xlabel('time  in  minutes') 
ylabel('grav  grad  z') 
grid 
pause 

if  printcon=l, 

print 

end 

% 

subplot(3 1  l),plot(time/60,aerox) 
title('aero  torques  vs  time') 
xlabel('time  in  minutes') 
ylabel('aero  x') 
grid 

subplot(3 1 2),plot(time/60,aeroy) 
xlabel('time  in  minutes') 
ylabel('aero  y') 
grid 

subplot(3 1 3),plot(time/60,aeroz) 
xlabel('time  in  minutes') 
ylabel('aero  z') 
grid 
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pause 


ifprintcon=l, 

print 

end 

end 

% 

clg 

subplot(4 1 1  ),plot(time/60 ,  wheelmomentum) 
title('momentum  wheel  characteristics  vs  time') 
xlabel('time  in  minutes') 
ylabel('wheel  momentum') 
grid 

subplot(4 1 2),plot(time/60,hdotwheel) 
xlabel('time  in  minutes') 
ylabel('h  dot  [wheel  torque]') 
grid 

subplot(4 1 3),plot(time/60,wheeldeltamomentum) 
xlabel('time  in  minutes') 
ylabel('delta  momentum') 
grid 

subplot(4 1 4),plot(time/60,  wheelspeed) 
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xlabel('time  in  minutes') 
ylabel('wheel  speed  in  rpm') 
grid 
pause 

if  printcon=l, 

print 

end 

% 

clg 

subplot(3 1  l),plot(time/60,magmomentsout(:,l)) 
title('magnetic  torque  rod  moments  vs  time') 
xlabel('time  in  minutes') 
ylabel('actual  mag  moment  x') 
grid 

subplot(3 1 2),plot(time/60,magmomentsout(:,2)) 
xlabel('time  in  minutes') 
ylabel('actual  mag  moment  y') 
grid 

subpIot(3 1 3),plot(time/60,magmomentsout(:,3)) 
xlabel('time  in  minutes') 
ylabel('actual  mag  moment  z') 
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grid 

pause 

if  printcon— 1, 
print 
end 
% 

subplot(3 1  l),plot(time/60,magtorquex) 
title('magnetic  torque  rod  torques  vs  time') 
xlabel('time  in  minutes') 
ylabel('mag  torque  x') 
grid 

subplot(3 1 2),plot(time/60,magtorquey) 
xlabel('time  in  minutes') 
ylabel('mag  torque  y') 
grid 

subplot(3 1 3),plot(time/60,magtorquez) 

xlabel('time  in  minutes') 

ylabel('mag  torque  z') 

grid 

pause 

if  printcon=l , 
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print 

end 

% 

subplot(3 1 1  ),plot(time/60,totaltorquein(:,  1 )) 
title('total  torques  each  axis  vs  time') 
xlabel('time  in  minutes') 
ylabel('total  torque  x  in  N-m') 
grid 

subplot(3 1 2),plot(time/60,totaltorquein(:,2)) 
xlabel('time  in  minutes') 
yIabel('total  torque  y  in  N-m') 
grid 

subplot(3 1 3),plot(time/60,totaltorquein(:,3)) 

xlabel('time  in  minutes') 

ylabel('total  torque  z  in  N-m') 

grid 

pause 

if  printcon=l, 
print 
end 
% 
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subplot(3 1  l),plot(time/60,disttorquex) 
title('dist  torques  each  axis  vs  time') 
xlabel('time  in  minutes') 
ylabel('dist  torque  x  in  N-m') 
grid 

subplot(3 1 2),plot(time/ 60,disttorquey) 
xlabel('time  in  minutes') 
ylabel('dist  torque  y  in  N-m') 
grid 

subplot(3 1 3),plot(time/60,disttorquez) 
xlabel('time  in  minutes') 
ylabel('dist  torque  z  in  N-m') 
grid 
pause 

if  printcon==l, 
print 
end 
% 

subplot(3 1  l),plot(time/60,bfldsensedbody(:,l)) 
title('sensed  b-field  vs  time') 
xlabel('time  in  minutes') 
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ylabel('sensed  b-field  x  dir') 
grid 

subplot(3 1 2),plot(time/60,bfldsensedbody(:,2)) 
xlabel('time  in  minutes') 
ylabel('sensed  b-field  y  dir') 
grid 

subplot(3 1 3),plot(time/60,bfldsensedbody(:,3)) 

xlabel('time  in  minutes') 

ylabel('sensed  b-field  z  dir') 

grid 

pause 

if  printcon==l, 
print 
end 
% 

%end 

% 

clg 

subplot(3 1  l),plot(time/60,omegaout(;,l)) 
title('body  rates') 
xlabel('time  in  minutes') 
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ylabel('omega  x') 
grid 

subplot(3 1 2),plot(time/60,omegaout(;,2)) 
xlabel('time  in  minutes') 
ylabel('omega  y') 
grid 

subplot(3 1 3),plot(time/60,omegaout(:,3)) 

xlabel('time  in  minutes') 

ylabel('omega  z') 

grid 

pause 

if  printcon=l, 

print 

end 

% 

clg 

subplot(3 1 1  ),plot(time/60,totalmomentum(:,  1 )) 

title('total  momentum') 

xlabel('time  in  minutes') 

ylabel('total  momentum  x') 

grid 
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subplot(3 1 2),plot(time/60,totalmomentum(:,2)) 
xlabel('time  in  minutes') 
ylabel('total  momentum  y') 
grid 

subplot(3 1 3),plot(time/60,totalmomentum(:,3)) 

xlabel('time  in  minutes') 

ylabel('totaI  momentum  z') 

grid 

pause 

if  printcon— 1 , 

print 

end 

% 

clg 

plot(time/60,pointerror*  1 80/pi, 'b') 

title('pointing  error  in  degrees  vs  time;  sun-oriented  control') 

xlabel('time  in  minutes') 

ylabel('pointing  error  to  sxm') 

grid 

hold 

plot(time/ 60,switchout,'b"') 
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pause 


if  printcon— 1, 

print 

end 
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